diff --git a/.github/workflows/check.yaml b/.github/workflows/check.yaml index abae2a1fb9..47c882c69b 100644 --- a/.github/workflows/check.yaml +++ b/.github/workflows/check.yaml @@ -27,7 +27,6 @@ jobs: secrets: REPO_GITHUB_TOKEN: ${{ secrets.REPO_GITHUB_TOKEN }} coverage: - if: github.event_name == 'pull_request' name: Coverage 📔 uses: insightsengineering/r.pkg.template/.github/workflows/test-coverage.yaml@main secrets: diff --git a/DESCRIPTION b/DESCRIPTION index 9454dd8849..98d6abfac1 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: tern Title: Create Common TLGs Used in Clinical Trials -Version: 0.7.8.9025 -Date: 2022-06-09 +Version: 0.7.9 +Date: 2022-08-16 Authors@R: c( person("Joe", "Zhu", , "joe.zhu@roche.com", role = c("aut", "cre")), person("Daniel", "Sabanés Bové", , "daniel.sabanes_bove@roche.com", role = "aut"), @@ -39,7 +39,7 @@ Imports: lifecycle, magrittr, methods, - nestcolor, + nestcolor (>= 0.0.1), rlang, scales, stats, @@ -51,8 +51,8 @@ Suggests: knitr, lattice, rmarkdown, - scda (>= 0.1.3), - scda.2022 (>= 0.1.1), + scda (>= 0.1.4), + scda.2022 (>= 0.1.2), stringr, testthat (>= 2.0) VignetteBuilder: @@ -96,6 +96,7 @@ Collate: 'd_pkparam.R' 'decorate_grob.R' 'deprec-nest_colors.R' + 'desctools_binom_diff.R' 'df_explicit_na.R' 'estimate_multinomial_rsp.R' 'estimate_proportion.R' diff --git a/NAMESPACE b/NAMESPACE index e90a06732c..8eece8045a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -59,6 +59,7 @@ export(as.rtable) export(color_palette) export(combine_groups) export(combine_levels) +export(combine_vectors) export(compare_vars) export(control_coxph) export(control_coxreg) @@ -111,6 +112,7 @@ export(extract_rsp_subgroups) export(extract_survival_biomarkers) export(extract_survival_subgroups) export(f_conf_level) +export(f_pval) export(fct_collapse_only) export(fct_discard) export(fct_explicit_na_if) @@ -252,6 +254,7 @@ export(score_occurrences_subtable) export(split_cols_by_groups) export(stack_grobs) export(stat_mean_ci) +export(stat_mean_pval) export(stat_median_ci) export(summarize_ancova) export(summarize_change) @@ -293,6 +296,7 @@ importFrom(methods,new) importFrom(rlang,"!!") importFrom(rlang,":=") importFrom(rlang,.data) +importFrom(stats,pchisq) importFrom(survival,Surv) importFrom(survival,coxph) importFrom(survival,strata) diff --git a/NEWS.md b/NEWS.md index 958ed04fb7..880a2b28a0 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,74 +1,54 @@ -# tern 0.7.8.9025 +# tern 0.7.9 ### Enhancements -* exported `draw_grob` function (needed by `enableRF`). -* `checkmate::assert()`, `checkmate::assert_true()`, and `checkmate::assert_false()` - use is kept to a minimum due to its ambiguous messages. -* `checkmate::assert_set_equal()` and `checkmate::assert_int()` have been changed - when checking for the length of vectors, either by checking the right input and - its length or by using `checkmate::assert_vector()`. -* Removed `assert_equal_length` as comparing lengths of multiple vectors can be - done without the need of a custom function. -* Implemented `nestcolor` in all examples with slight refactoring to `g_km`, `g_ipp`, +* Added `DescTools` `BinomDiffCI` function within `tern`. +* Exported `draw_grob` function. +* Added new parameter to `summarize_logistic` to specify which pivoted value should be removed during analysis. +* Updated `s_coxph_pairwise` to generate log-rank p-value using original log-rank test instead + of Cox Proportional-Hazards Model. +* Added more verbose warnings from `as_factor_keep_attributes` for better `test_examples()` readability. +* Implemented `nestcolor` in all examples, with slight refactoring to `g_km`, `g_ipp`, `g_waterfall`, `g_step`, `g_lineplot`, and `g_forest`. +* Added `stat_mean_pval` function to calculate the p-value of the mean as a new summary statistic. +* New statistic `mean_se` (mean with standard error) for `summarize_variables()` + and related functions. ### Migration from `assertthat` to `checkmate` -* complete substitution of `assertthat` calls with `checkmate`. -* `assert_df_with_factors`, `assert_equal_length`, and `assert_proportion_value` - made internals. -* removed `has_tabletree_colnames`. -* removed `is_quantiles_vector` and `all_elements_in_ref` (substituted by - `checkmate::subset()`). -* removed `is_proportion_vector` and adapted `is_proportion` to `checkmate` style - with `assert_proportion_value`. -* `is_equal_length` is now `assert_equal_length`. -* removed `is_df_with_no_na_level` by adding parameters to `assert_df_with_variables`. -* removed `is_df_with_nlevels_factor` by adding parameters to `assert_df_with_factors`. -* `is_character_or_factor` replaced by `checkmate::assert_multi_class()`. -* `is_nonnegative_count` removed and substituted by `checkmate::assert_count()`. -* `assert_list_of_variables` instead of `is_variables` (made internal). -* `assert_df_with_variables` instead of `is_df_with_variables` (made internal). -* `assert_df_with_factors` instead of `is_df_with_factors` (made internal). -* `assert_valid_factor` instead of `is_valid_factor`. -* removed `is_valid_character`. -* renamed `assertthat.R` into a more appropriate `utils_checkmate.R`. -* renamed `test-assertthat.R` into `test-utils_checkmate.R`. +* Substituted all `assertthat` calls with `checkmate`. +* Implemented `checkmate::assert_vector()`, `checkmate::assert_set_equal()`, and + `checkmate::assert_int()` to validate vector type, length, and input. +* Made `assert_df_with_factors` and `assert_proportion_value` internal functions. +* Removed assertion functions `all_elements_in_ref`, `is_df_with_nlevels_factor`, + `is_df_with_no_na_level`, `is_proportion`, `is_proportion_vector`, `is_quantiles_vector`, + `is_valid_character`, + `assert_character_or_factor`, `assert_equal_length` and `has_tabletree_colnames`. +* Renamed `is_proportion` to `assert_proportion_value`. +* Renamed `is_equal_length` to `assert_equal_length`. +* Refactored `assert_df_with_variables` to replace `is_df_with_no_na_level`. +* Refactored `assert_df_with_factors` to replace `is_df_with_nlevels_factor`. +* Replaced usage of `is_quantiles_vector` and `all_elements_in_ref` with `checkmate::subset()`. +* Replaced `is_character_or_factor` with `checkmate::assert_multi_class()`. +* Replaced `is_nonnegative_count` with `checkmate::assert_count()`. +* Replaced `is_variables` with `assert_list_of_variables`. +* Replaced `is_df_with_variables` with `assert_df_with_variables`. +* Replaced `is_df_with_factors` with `assert_df_with_factors`. +* Replaced `is_valid_factor` with `assert_valid_factor`. +* Renamed `assertthat.R` to `utils_checkmate.R`. +* Renamed `test-assertthat.R` to `test-utils_checkmate.R`. -### Fix -* Fixing error coming from comparing factors vector to characters vector. -* `fct_collapse_only`, `fct_collapse_only`, and `month2day`/`day2month` reverted - to export. -* Fixed test for `cut_quantile_bins` with empty vector. This is not a good standard - input. -* Warnings from `as_factor_keep_attributes` are now in verbose for better - `test_examples()` readability -* Renaming `rtables.R` as confusing file name due to the package dependence. -* Renamed files to respect the main documented function and fixed file in `R/` - folder that does not comply with the standard (starting with the word test). -* Extracted `cox_regression_inter` from `cox_regression`. -* `d_` and `h_` functions are all reverted or confirmed as export functions to - allow future users to utilize their added flexibility. -* Fixing bug related to error flag for empty strings coming from `rtables` split +### Bug Fixes +* Identified bug in `prop_diff` functions. Coding of responses have been corrected (TRUE is a success). +* Fixed error coming from comparing factors vector to characters vector. +* Fixed empty vector test for `cut_quantile_bins`. +* Fixed bug related to error flag for empty strings coming from `rtables` split functions. Creation of `replace_emptys_with_na` to replace empty strings with custom strings across data.frame (this can be merged with `df_explicit_na`). -* Add of new parameter for `summarize_logistic` that specify which pivoted value - is meant to be removed during the analysis. -* Renamed `estimate_incidence_rate.R` into `incidence_rate.R` to match the - documentation grouping name. -* Extracted `control_incidence_rate` into a separated file as it was exported and - documented separately. -* Added `@md` and removed `@order` from `incidence_rate.R`. The presented order - is automatically the final order in the documentation. There is no specific - need to add this tag. Modified examples accordingly. -* Fixed warnings occurring in example tests -* Removed `tern:::` prefix and added `dontrun` to internal function examples. -* Enhanced `s_coxph_pairwise` with generating log-rank p value by original - log-rank test, instead of Cox Proportional-Hazards Model. -* Fixed bug in `s_ancova`, avoiding error when the first level of the arm factor - is not the control arm. +* Fixed warnings occurring in example tests. +* Fixed internal function examples errors by removing `tern:::` prefix and added `dontrun` to internal function examples. +* Fixed bug in `s_ancova` causing an error when the first level of the arm factor is not the control arm. -### Documentation and NAMESPACE polishing -* Added stable badge for: +### Documentation and NAMESPACE Polishing +* Added stable badges for: - `count_abnormal_by_marked` (reference to `abnormal_by_marked`), `count_abnormal_lab_worsen_by_baseline` and `h_adlb_worsen` (reference to `abnormal_by_worst_grade_worsen_from_baseline`), `count_abnormal_by_worst_grade` @@ -81,7 +61,7 @@ `a_summary` (reference to `summarize_variables`) and kept the S3 method tree. - `summarize_patients_exposure_in_cols`, `summarize_num_patients` with `s_num_patients`, `s_num_patients_content`, `summarize_num_patients`. - - `count_cumulative`, `count_missed_doses`, `count_patients_events_in_cols`, `summarize_colvars`, `summarize_change`, `summarize_ancova`,`as.rtable`, `color_palette`, `add_footnotes` (note, this function is defined in two different files: footnotes and g_forest). + - `count_cumulative`, `count_missed_doses`, `count_patients_events_in_cols`, `summarize_colvars`, `summarize_change`, `summarize_ancova`,`as.rtable`, `color_palette`, `add_footnotes`. - (statistical function controls) `control_coxreg`, `control_coxph`, `control_incidence_rate`, `control_lineplot_vars`, `control_surv_time`, `control_surv_timepoint`, `control_logisitic`, `control_step`. @@ -89,16 +69,15 @@ `tabulate_rsp_biomarkers`, `keep_rows`, `keep_content_rows`, `has_count_in_any_col`, `has_fraction_in_cols`, `has_fraction_in_any_col`, `has_fractions_difference`, `test_proportion_diff`, `pairwise`, `logistic_regression`, - `estimate_incidence_rate`, `control_incidence_rate` (another file), reference to - `incidence_rate`, `cut_quantile_bins`, `d_pkparam` (note: d function but required by tern_ex file), + `estimate_incidence_rate`, `control_incidence_rate` (reference to + `incidence_rate`), `cut_quantile_bins`, `d_pkparam`, `estimate_multinomial_rsp`, `decorate_grob_set`, `extreme_format`, `fit_rsp_step`, - `fit_survival_step`, `footnotes`, `footnotes-set` (note: this function is defined twice in the footnotes file), + `fit_survival_step`, `footnotes`, `footnotes-set`, `format_count_fraction`, `format_fraction_threshold`, `formatting_functions`, `format_fraction`, `combination_function` (S4 method), `compare_variables` (S3 method), `h_stack_by_baskets`, `h_pkparam_sort`, `h_adsl_adlb_merge_using_worst_flag`, `kaplan_meier`. - -* Internal keywords added, export removed, `_pkgdown.yml` polished and `tern:::` for - tests, examples, and vignettes when present for the following functions: +* Internal keywords added, export removed, `_pkgdown.yml` updated, and `tern:::` added for + tests/examples/vignettes where present for the following functions: - (chain functions, reference to `abnormal_by_marked`) `s_count_abnormal_by_marked`, `a_count_abnormal_by_marked`. - (chain functions, reference to `abnormal_by_worst_grade_worsen_from_baseline`) @@ -133,18 +112,26 @@ `h_content_first_row`, `a_response_subgroups`, `range_noinf`, `has_count_in_cols`, `has_counts_difference`, `prop_chisq`, `prop_cmh`, `prop_schouten`, `prop_fisher`, `s_test_proportion_diff`, `a_test_proportion_diff`,`h_split_param`, - `fct_collapse_only`, `fct_discard`, `fct_explicit_na_if`. - -* Deprecated badge added to `g_mmrm` -* Removed `tern:::` prefix from internal function uses in tests + `fct_discard`, `fct_explicit_na_if`. +* Added deprecated badge to `g_mmrm`. +* Removed `tern:::` prefix from internal function uses in tests. +* Reverted `fct_collapse_only`, `month2day`, and `day2month` to export. +* Reverted all `d_` and `h_` functions to export, allowing users to utilize their flexibility. +* Renamed `rtables.R` as file name is confusing due to the package dependence. +* Renamed files to reflect main documented function and fixed `R` files starting with the word test + not complying with the standard. +* Extracted `cox_regression_inter` into a separate file from `cox_regression`. +* Renamed `estimate_incidence_rate.R` to `incidence_rate.R` to match the documentation grouping name. +* Extracted `control_incidence_rate` into a separate file as it was exported and documented separately. +* Added `@md` and removed `@order` from `incidence_rate.R`. Modified examples accordingly. +* Removed hyperlink from `prop_schouten` function documentation. ### Miscellaneous * Deprecated the `color_palette` function with `nestcolor::color_palette` and removed `color_palette_core`. * Deprecated `h_set_nest_theme()` in favor of `nestcolor::theme_nest()`. * Removed deprecated `mmrm` functions: `fit_mmrm`, `g_mmrm_diagnostic`, `g_mmrm_lsmeans`, `as.rtable.mmrm`, `h_mmrm_fixed`, `h_mmrm_cov`, `h_mmrm_diagnostic`, `tidy.mmrm`, `s_mmrm_lsmeans`, `s_mmrm_lsmeans_single`, `summarize_lsmeans`. -* Changed function names `arm` to `study_arm` and `extract` to `extract_by_name`. - +* Renamed functions `arm` to `study_arm` and `extract` to `extract_by_name`. # tern 0.7.8 @@ -167,7 +154,7 @@ tests, examples, and vignettes when present for the following functions: - (helper functions) `h_format_row`, `h_map_for_count_abnormal` - (utils functions) `make_names`, `month2day`, `day2month` - `empty_vector_if_na`, `combine_vectors`, `aesi_label`, + `empty_vector_if_na`, `aesi_label`, `n_available`, `format_xx`, `arm`. - `count_values_funs`, `prop_difference`, `combine_counts`. - (chain functions) `s_count_abnormal`, `a_count_abnormal`. diff --git a/R/desctools_binom_diff.R b/R/desctools_binom_diff.R new file mode 100644 index 0000000000..26fa77601d --- /dev/null +++ b/R/desctools_binom_diff.R @@ -0,0 +1,747 @@ +#' Confidence Intervals for a Difference of Binomials +#' +#' @description `r lifecycle::badge("experimental")` +#' +#' Several confidence intervals for the difference between proportions. +#' +#' @param grp (`factor`)\cr +#' vector assigning observations to one out of two groups +#' (e.g. reference and treatment group). +#' @name desctools_binom + +NULL + +#' Recycle list of parameters +#' +#' @describeIn desctools_binom This function recycles all supplied elements to the maximal dimension. +#' +#' @keywords internal +h_recycle <- function(...) { + lst <- list(...) + maxdim <- max(lengths(lst)) + res <- lapply(lst, rep, length.out = maxdim) + attr(res, "maxdim") <- maxdim + return(res) +} + +#' @describeIn desctools_binom Several Confidence Intervals for the difference between proportions. +#' +#' +#' @return A named list of 3 values: +#' - `est`: estimate of proportion difference. +#' - `lwrci`: estimate of lower end of the confidence interval +#' - `upci`: estimate of upper end of the confidence interval. +#' @examples +#' # Internal function - +#' \dontrun{ +#' +#' set.seed(2) +#' rsp <- sample(c(TRUE, FALSE), replace = TRUE, size = 20) +#' grp <- factor(c(rep("A", 10), rep("B", 10))) +#' tbl <- table(grp, factor(rsp, levels = c(TRUE, FALSE))) +#' desctools_binom(tbl[1], sum(tbl[1], tbl[3]), tbl[2], sum(tbl[2], tbl[4]), conf.level = 0.90, method = "waldcc") +#' } +#' +#' @keywords internal +desctools_binom <- function(x1, n1, x2, n2, conf.level = 0.95, sides = c( + "two.sided", + "left", "right" + ), method = c( + "ac", "wald", "waldcc", "score", + "scorecc", "mn", "mee", "blj", "ha", "hal", "jp" + )) { + if (missing(sides)) { + sides <- match.arg(sides) + } + if (missing(method)) { + method <- match.arg(method) + } + iBinomDiffCI <- function(x1, n1, x2, n2, conf.level, sides, + method) { + if (sides != "two.sided") { + conf.level <- 1 - 2 * (1 - conf.level) + } + alpha <- 1 - conf.level + kappa <- stats::qnorm(1 - alpha / 2) + p1.hat <- x1 / n1 + p2.hat <- x2 / n2 + est <- p1.hat - p2.hat + switch(method, + wald = { + vd <- p1.hat * (1 - p1.hat) / n1 + p2.hat * (1 - p2.hat) / n2 + term2 <- kappa * sqrt(vd) + CI.lower <- max(-1, est - term2) + CI.upper <- min(1, est + term2) + }, + waldcc = { + vd <- p1.hat * (1 - p1.hat) / n1 + p2.hat * (1 - p2.hat) / n2 + term2 <- kappa * sqrt(vd) + term2 <- term2 + 0.5 * (1 / n1 + 1 / n2) + CI.lower <- max(-1, est - term2) + CI.upper <- min(1, est + term2) + }, + ac = { + n1 <- n1 + 2 + n2 <- n2 + 2 + x1 <- x1 + 1 + x2 <- x2 + 1 + p1.hat <- x1 / n1 + p2.hat <- x2 / n2 + est1 <- p1.hat - p2.hat + vd <- p1.hat * (1 - p1.hat) / n1 + p2.hat * (1 - p2.hat) / n2 + term2 <- kappa * sqrt(vd) + CI.lower <- max(-1, est1 - term2) + CI.upper <- min(1, est1 + term2) + }, + exact = { + CI.lower <- NA + CI.upper <- NA + }, + score = { + w1 <- desctools_binomci( + x = x1, n = n1, conf.level = conf.level, + method = "wilson" + ) + w2 <- desctools_binomci( + x = x2, n = n2, conf.level = conf.level, + method = "wilson" + ) + l1 <- w1[2] + u1 <- w1[3] + l2 <- w2[2] + u2 <- w2[3] + CI.lower <- est - kappa * sqrt(l1 * (1 - l1) / n1 + + u2 * (1 - u2) / n2) + CI.upper <- est + kappa * sqrt(u1 * (1 - u1) / n1 + + l2 * (1 - l2) / n2) + }, + scorecc = { + w1 <- desctools_binomci( + x = x1, n = n1, conf.level = conf.level, + method = "wilsoncc" + ) + w2 <- desctools_binomci( + x = x2, n = n2, conf.level = conf.level, + method = "wilsoncc" + ) + l1 <- w1[2] + u1 <- w1[3] + l2 <- w2[2] + u2 <- w2[3] + CI.lower <- max(-1, est - sqrt((p1.hat - l1)^2 + + (u2 - p2.hat)^2)) + CI.upper <- min(1, est + sqrt((u1 - p1.hat)^2 + (p2.hat - + l2)^2)) + }, + mee = { + .score <- function(p1, n1, p2, n2, dif) { + if (dif > 1) dif <- 1 + if (dif < -1) dif <- -1 + diff <- p1 - p2 - dif + if (abs(diff) == 0) { + res <- 0 + } else { + t <- n2 / n1 + a <- 1 + t + b <- -(1 + t + p1 + t * p2 + dif * (t + 2)) + c <- dif * dif + dif * (2 * p1 + t + 1) + p1 + + t * p2 + d <- -p1 * dif * (1 + dif) + v <- (b / a / 3)^3 - b * c / (6 * a * a) + d / a / 2 + if (abs(v) < .Machine$double.eps) v <- 0 + s <- sqrt((b / a / 3)^2 - c / a / 3) + u <- ifelse(v > 0, 1, -1) * s + w <- (3.141592654 + acos(v / u^3)) / 3 + p1d <- 2 * u * cos(w) - b / a / 3 + p2d <- p1d - dif + n <- n1 + n2 + res <- (p1d * (1 - p1d) / n1 + p2d * (1 - p2d) / n2) + } + return(sqrt(res)) + } + pval <- function(delta) { + z <- (est - delta) / .score( + p1.hat, n1, p2.hat, + n2, delta + ) + 2 * min(stats::pnorm(z), 1 - stats::pnorm(z)) + } + CI.lower <- max(-1, stats::uniroot(function(delta) { + pval(delta) - + alpha + }, interval = c(-1 + 1e-06, est - 1e-06))$root) + CI.upper <- min(1, stats::uniroot(function(delta) { + pval(delta) - + alpha + }, interval = c(est + 1e-06, 1 - 1e-06))$root) + }, + blj = { + p1.dash <- (x1 + 0.5) / (n1 + 1) + p2.dash <- (x2 + 0.5) / (n2 + 1) + vd <- p1.dash * (1 - p1.dash) / n1 + p2.dash * (1 - + p2.dash) / n2 + term2 <- kappa * sqrt(vd) + est.dash <- p1.dash - p2.dash + CI.lower <- max(-1, est.dash - term2) + CI.upper <- min(1, est.dash + term2) + }, + ha = { + term2 <- 1 / (2 * min(n1, n2)) + kappa * sqrt(p1.hat * + (1 - p1.hat) / (n1 - 1) + p2.hat * (1 - p2.hat) / (n2 - + 1)) + CI.lower <- max(-1, est - term2) + CI.upper <- min(1, est + term2) + }, + mn = { + .conf <- function(x1, n1, x2, n2, z, lower = FALSE) { + p1 <- x1 / n1 + p2 <- x2 / n2 + p.hat <- p1 - p2 + dp <- 1 + ifelse(lower, 1, -1) * p.hat + i <- 1 + while (i <= 50) { + dp <- 0.5 * dp + y <- p.hat + ifelse(lower, -1, 1) * dp + score <- .score(p1, n1, p2, n2, y) + if (score < z) { + p.hat <- y + } + if ((dp < 1e-07) || (abs(z - score) < 1e-06)) { + (break)() + } else { + i <- i + + 1 + } + } + return(y) + } + .score <- function(p1, n1, p2, n2, dif) { + diff <- p1 - p2 - dif + if (abs(diff) == 0) { + res <- 0 + } else { + t <- n2 / n1 + a <- 1 + t + b <- -(1 + t + p1 + t * p2 + dif * (t + 2)) + c <- dif * dif + dif * (2 * p1 + t + 1) + p1 + + t * p2 + d <- -p1 * dif * (1 + dif) + v <- (b / a / 3)^3 - b * c / (6 * a * a) + d / a / 2 + s <- sqrt((b / a / 3)^2 - c / a / 3) + u <- ifelse(v > 0, 1, -1) * s + w <- (3.141592654 + acos(v / u^3)) / 3 + p1d <- 2 * u * cos(w) - b / a / 3 + p2d <- p1d - dif + n <- n1 + n2 + var <- (p1d * (1 - p1d) / n1 + p2d * (1 - p2d) / n2) * + n / (n - 1) + res <- diff^2 / var + } + return(res) + } + z <- stats::qchisq(conf.level, 1) + CI.lower <- max(-1, .conf(x1, n1, x2, n2, z, TRUE)) + CI.upper <- min(1, .conf(x1, n1, x2, n2, z, FALSE)) + }, + beal = { + a <- p1.hat + p2.hat + b <- p1.hat - p2.hat + u <- ((1 / n1) + (1 / n2)) / 4 + v <- ((1 / n1) - (1 / n2)) / 4 + V <- u * ((2 - a) * a - b^2) + 2 * v * (1 - a) * + b + z <- stats::qchisq(p = 1 - alpha / 2, df = 1) + A <- sqrt(z * (V + z * u^2 * (2 - a) * a + z * v^2 * + (1 - a)^2)) + B <- (b + z * v * (1 - a)) / (1 + z * u) + CI.lower <- max(-1, B - A / (1 + z * u)) + CI.upper <- min(1, B + A / (1 + z * u)) + }, + hal = { + psi <- (p1.hat + p2.hat) / 2 + u <- (1 / n1 + 1 / n2) / 4 + v <- (1 / n1 - 1 / n2) / 4 + z <- kappa + theta <- ((p1.hat - p2.hat) + z^2 * v * (1 - 2 * + psi)) / (1 + z^2 * u) + w <- z / (1 + z^2 * u) * sqrt(u * (4 * psi * (1 - psi) - + (p1.hat - p2.hat)^2) + 2 * v * (1 - 2 * psi) * + (p1.hat - p2.hat) + 4 * z^2 * u^2 * (1 - psi) * + psi + z^2 * v^2 * (1 - 2 * psi)^2) + c(theta + w, theta - w) + CI.lower <- max(-1, theta - w) + CI.upper <- min(1, theta + w) + }, + jp = { + psi <- 0.5 * ((x1 + 0.5) / (n1 + 1) + (x2 + 0.5) / (n2 + + 1)) + u <- (1 / n1 + 1 / n2) / 4 + v <- (1 / n1 - 1 / n2) / 4 + z <- kappa + theta <- ((p1.hat - p2.hat) + z^2 * v * (1 - 2 * + psi)) / (1 + z^2 * u) + w <- z / (1 + z^2 * u) * sqrt(u * (4 * psi * (1 - psi) - + (p1.hat - p2.hat)^2) + 2 * v * (1 - 2 * psi) * + (p1.hat - p2.hat) + 4 * z^2 * u^2 * (1 - psi) * + psi + z^2 * v^2 * (1 - 2 * psi)^2) + c(theta + w, theta - w) + CI.lower <- max(-1, theta - w) + CI.upper <- min(1, theta + w) + }, + ) + ci <- c( + est = est, lwr.ci = min(CI.lower, CI.upper), + upr.ci = max(CI.lower, CI.upper) + ) + if (sides == "left") { + ci[3] <- 1 + } else if (sides == "right") { + ci[2] <- -1 + } + return(ci) + } + method <- match.arg(arg = method, several.ok = TRUE) + sides <- match.arg(arg = sides, several.ok = TRUE) + lst <- h_recycle( + x1 = x1, n1 = n1, x2 = x2, n2 = n2, conf.level = conf.level, + sides = sides, method = method + ) + res <- t(sapply(1:attr(lst, "maxdim"), function(i) { + iBinomDiffCI( + x1 = lst$x1[i], + n1 = lst$n1[i], x2 = lst$x2[i], n2 = lst$n2[i], conf.level = lst$conf.level[i], + sides = lst$sides[i], method = lst$method[i] + ) + })) + lgn <- h_recycle(x1 = if (is.null(names(x1))) { + paste("x1", seq_along(x1), sep = ".") + } else { + names(x1) + }, n1 = if (is.null(names(n1))) { + paste("n1", seq_along(n1), sep = ".") + } else { + names(n1) + }, x2 = if (is.null(names(x2))) { + paste("x2", seq_along(x2), sep = ".") + } else { + names(x2) + }, n2 = if (is.null(names(n2))) { + paste("n2", seq_along(n2), sep = ".") + } else { + names(n2) + }, conf.level = conf.level, sides = sides, method = method) + xn <- apply(as.data.frame(lgn[sapply(lgn, function(x) { + length(unique(x)) != + 1 + })]), 1, paste, collapse = ":") + rownames(res) <- xn + return(res) +} + +#' @describeIn desctools_binom Compute confidence intervals for binomial proportions. +#' +#' @param x \cr number of successes +#' @param n \cr number of trials +#' @param conf.level \cr confidence level, defaults to 0.95 +#' @param sides \cr a character string specifying the side of the confidence interval. Must be one of "two-sided" (default), "left" or "right". +#' @param method \cr character string specifying which method to use. Can be one out of: "wald", "wilson", "wilsoncc", "agresti-coull", "jeffreys", "modified wilson", "modified jeffreys", "clopper-pearson", "arcsine. +#' ", "logit", "witting", "pratt", "midp", "lik" and "blaker" +#' +#' +#' @return A matric with 3 columns containing: +#' - `est`: estimate of proportion difference. +#' - `lwrci`: lower end of the confidence interval +#' - `upci`: upper end of the confidence interval. +#' +#' @keywords internal +desctools_binomci <- function(x, n, conf.level = 0.95, sides = c( + "two.sided", "left", + "right" + ), method = c( + "wilson", "wald", "waldcc", "agresti-coull", + "jeffreys", "modified wilson", "wilsoncc", "modified jeffreys", + "clopper-pearson", "arcsine", "logit", "witting", "pratt", + "midp", "lik", "blaker" + ), rand = 123, tol = 1e-05) { + if (missing(method)) { + method <- "wilson" + } + if (missing(sides)) { + sides <- "two.sided" + } + iBinomCI <- function(x, n, conf.level = 0.95, sides = c( + "two.sided", + "left", "right" + ), method = c( + "wilson", "wilsoncc", "wald", + "waldcc", "agresti-coull", "jeffreys", "modified wilson", + "modified jeffreys", "clopper-pearson", "arcsine", "logit", + "witting", "pratt", "midp", "lik", "blaker" + ), rand = 123, + tol = 1e-05) { + if (length(x) != 1) { + stop("'x' has to be of length 1 (number of successes)") + } + if (length(n) != 1) { + stop("'n' has to be of length 1 (number of trials)") + } + if (length(conf.level) != 1) { + stop("'conf.level' has to be of length 1 (confidence level)") + } + if (conf.level < 0.5 | conf.level > 1) { + stop("'conf.level' has to be in [0.5, 1]") + } + sides <- match.arg(sides, choices = c( + "two.sided", "left", + "right" + ), several.ok = FALSE) + if (sides != "two.sided") { + conf.level <- 1 - 2 * (1 - conf.level) + } + alpha <- 1 - conf.level + kappa <- stats::qnorm(1 - alpha / 2) + p.hat <- x / n + q.hat <- 1 - p.hat + est <- p.hat + switch(match.arg(arg = method, choices = c( + "wilson", + "wald", "waldcc", "wilsoncc", "agresti-coull", "jeffreys", + "modified wilson", "modified jeffreys", "clopper-pearson", + "arcsine", "logit", "witting", "pratt", "midp", "lik", + "blaker" + )), + wald = { + term2 <- kappa * sqrt(p.hat * q.hat) / sqrt(n) + CI.lower <- max(0, p.hat - term2) + CI.upper <- min(1, p.hat + term2) + }, + waldcc = { + term2 <- kappa * sqrt(p.hat * q.hat) / sqrt(n) + term2 <- term2 + 1 / (2 * n) + CI.lower <- max(0, p.hat - term2) + CI.upper <- min(1, p.hat + term2) + }, + wilson = { + term1 <- (x + kappa^2 / 2) / (n + kappa^2) + term2 <- kappa * sqrt(n) / (n + kappa^2) * sqrt(p.hat * + q.hat + kappa^2 / (4 * n)) + CI.lower <- max(0, term1 - term2) + CI.upper <- min(1, term1 + term2) + }, + wilsoncc = { + lci <- (2 * x + kappa^2 - 1 - kappa * sqrt(kappa^2 - + 2 - 1 / n + 4 * p.hat * (n * q.hat + 1))) / (2 * + (n + kappa^2)) + uci <- (2 * x + kappa^2 + 1 + kappa * sqrt(kappa^2 + + 2 - 1 / n + 4 * p.hat * (n * q.hat - 1))) / (2 * + (n + kappa^2)) + CI.lower <- max(0, ifelse(p.hat == 0, 0, lci)) + CI.upper <- min(1, ifelse(p.hat == 1, 1, uci)) + }, + `agresti-coull` = { + x.tilde <- x + kappa^2 / 2 + n.tilde <- n + kappa^2 + p.tilde <- x.tilde / n.tilde + q.tilde <- 1 - p.tilde + est <- p.tilde + term2 <- kappa * sqrt(p.tilde * q.tilde) / sqrt(n.tilde) + CI.lower <- max(0, p.tilde - term2) + CI.upper <- min(1, p.tilde + term2) + }, + jeffreys = { + if (x == 0) { + CI.lower <- 0 + } else { + CI.lower <- stats::qbeta( + alpha / 2, + x + 0.5, n - x + 0.5 + ) + } + if (x == n) { + CI.upper <- 1 + } else { + CI.upper <- stats::qbeta(1 - + alpha / 2, x + 0.5, n - x + 0.5) + } + }, + `modified wilson` = { + term1 <- (x + kappa^2 / 2) / (n + kappa^2) + term2 <- kappa * sqrt(n) / (n + kappa^2) * sqrt(p.hat * + q.hat + kappa^2 / (4 * n)) + if ((n <= 50 & x %in% c(1, 2)) | (n >= 51 & x %in% + c(1:3))) { + CI.lower <- 0.5 * stats::qchisq(alpha, 2 * + x) / n + } else { + CI.lower <- max(0, term1 - term2) + } + if ((n <= 50 & x %in% c(n - 1, n - 2)) | (n >= 51 & + x %in% c(n - (1:3)))) { + CI.upper <- 1 - 0.5 * stats::qchisq( + alpha, + 2 * (n - x) + ) / n + } else { + CI.upper <- min(1, term1 + + term2) + } + }, + `modified jeffreys` = { + if (x == n) { + CI.lower <- (alpha / 2)^(1 / n) + } else { + if (x <= 1) { + CI.lower <- 0 + } else { + CI.lower <- stats::qbeta( + alpha / 2, + x + 0.5, n - x + 0.5 + ) + } + } + if (x == 0) { + CI.upper <- 1 - (alpha / 2)^(1 / n) + } else { + if (x >= n - 1) { + CI.upper <- 1 + } else { + CI.upper <- stats::qbeta(1 - + alpha / 2, x + 0.5, n - x + 0.5) + } + } + }, + `clopper-pearson` = { + CI.lower <- stats::qbeta(alpha / 2, x, n - x + 1) + CI.upper <- stats::qbeta(1 - alpha / 2, x + 1, n - x) + }, + arcsine = { + p.tilde <- (x + 0.375) / (n + 0.75) + est <- p.tilde + CI.lower <- sin(asin(sqrt(p.tilde)) - 0.5 * kappa / sqrt(n))^2 + CI.upper <- sin(asin(sqrt(p.tilde)) + 0.5 * kappa / sqrt(n))^2 + }, + logit = { + lambda.hat <- log(x / (n - x)) + V.hat <- n / (x * (n - x)) + lambda.lower <- lambda.hat - kappa * sqrt(V.hat) + lambda.upper <- lambda.hat + kappa * sqrt(V.hat) + CI.lower <- exp(lambda.lower) / (1 + exp(lambda.lower)) + CI.upper <- exp(lambda.upper) / (1 + exp(lambda.upper)) + }, + witting = { + set.seed(rand) + x.tilde <- x + stats::runif(1, min = 0, max = 1) + pbinom.abscont <- function(q, size, prob) { + v <- trunc(q) + term1 <- stats::pbinom(v - 1, size = size, prob = prob) + term2 <- (q - v) * stats::dbinom(v, size = size, prob = prob) + return(term1 + term2) + } + qbinom.abscont <- function(p, size, x) { + fun <- function(prob, size, x, p) { + pbinom.abscont(x, size, prob) - p + } + stats::uniroot(fun, + interval = c(0, 1), size = size, + x = x, p = p + )$root + } + CI.lower <- qbinom.abscont(1 - alpha, size = n, x = x.tilde) + CI.upper <- qbinom.abscont(alpha, size = n, x = x.tilde) + }, + pratt = { + if (x == 0) { + CI.lower <- 0 + CI.upper <- 1 - alpha^(1 / n) + } else if (x == 1) { + CI.lower <- 1 - (1 - alpha / 2)^(1 / n) + CI.upper <- 1 - (alpha / 2)^(1 / n) + } else if (x == (n - 1)) { + CI.lower <- (alpha / 2)^(1 / n) + CI.upper <- (1 - alpha / 2)^(1 / n) + } else if (x == n) { + CI.lower <- alpha^(1 / n) + CI.upper <- 1 + } else { + z <- stats::qnorm(1 - alpha / 2) + A <- ((x + 1) / (n - x))^2 + B <- 81 * (x + 1) * (n - x) - 9 * n - 8 + C <- (0 - 3) * z * sqrt(9 * (x + 1) * (n - x) * + (9 * n + 5 - z^2) + n + 1) + D <- 81 * (x + 1)^2 - 9 * (x + 1) * (2 + z^2) + + 1 + E <- 1 + A * ((B + C) / D)^3 + CI.upper <- 1 / E + A <- (x / (n - x - 1))^2 + B <- 81 * x * (n - x - 1) - 9 * n - 8 + C <- 3 * z * sqrt(9 * x * (n - x - 1) * (9 * + n + 5 - z^2) + n + 1) + D <- 81 * x^2 - 9 * x * (2 + z^2) + 1 + E <- 1 + A * ((B + C) / D)^3 + CI.lower <- 1 / E + } + }, + midp = { + f.low <- function(pi, x, n) { + 1 / 2 * stats::dbinom(x, size = n, prob = pi) + stats::pbinom(x, + size = n, prob = pi, lower.tail = FALSE + ) - + (1 - conf.level) / 2 + } + f.up <- function(pi, x, n) { + 1 / 2 * stats::dbinom(x, size = n, prob = pi) + stats::pbinom(x - + 1, size = n, prob = pi) - (1 - conf.level) / 2 + } + CI.lower <- 0 + CI.upper <- 1 + if (x != 0) { + CI.lower <- stats::uniroot(f.low, + interval = c(0, p.hat), + x = x, n = n + )$root + } + if (x != n) { + CI.upper <- stats::uniroot(f.up, interval = c( + p.hat, + 1 + ), x = x, n = n)$root + } + }, + lik = { + CI.lower <- 0 + CI.upper <- 1 + z <- stats::qnorm(1 - alpha * 0.5) + tol <- .Machine$double.eps^0.5 + BinDev <- function(y, x, mu, wt, bound = 0, tol = .Machine$double.eps^0.5, + ...) { + ll.y <- ifelse(y %in% c(0, 1), 0, stats::dbinom(x, wt, + y, + log = TRUE + )) + ll.mu <- ifelse(mu %in% c(0, 1), 0, stats::dbinom(x, + wt, mu, + log = TRUE + )) + res <- ifelse(abs(y - mu) < tol, 0, sign(y - + mu) * sqrt(-2 * (ll.y - ll.mu))) + return(res - bound) + } + if (x != 0 && tol < p.hat) { + CI.lower <- if (BinDev( + tol, x, p.hat, n, -z, + tol + ) <= 0) { + stats::uniroot( + f = BinDev, interval = c(tol, if (p.hat < + tol || p.hat == 1) { + 1 - tol + } else { + p.hat + }), bound = -z, + x = x, mu = p.hat, wt = n + )$root + } + } + if (x != n && p.hat < (1 - tol)) { + CI.upper <- if (BinDev(y = 1 - tol, x = x, mu = ifelse(p.hat > + 1 - tol, tol, p.hat), wt = n, bound = z, tol = tol) < + 0) { + CI.lower <- if (BinDev( + tol, x, if (p.hat < + tol || p.hat == 1) { + 1 - tol + } else { + p.hat + }, n, + -z, tol + ) <= 0) { + stats::uniroot( + f = BinDev, interval = c(tol, p.hat), + bound = -z, x = x, mu = p.hat, wt = n + )$root + } + } else { + stats::uniroot( + f = BinDev, interval = c(if (p.hat > + 1 - tol) { + tol + } else { + p.hat + }, 1 - tol), bound = z, + x = x, mu = p.hat, wt = n + )$root + } + } + }, + blaker = { + acceptbin <- function(x, n, p) { + p1 <- 1 - stats::pbinom(x - 1, n, p) + p2 <- stats::pbinom(x, n, p) + a1 <- p1 + stats::pbinom(stats::qbinom(p1, n, p) - 1, n, p) + a2 <- p2 + 1 - stats::pbinom( + stats::qbinom(1 - p2, n, p), n, + p + ) + return(min(a1, a2)) + } + CI.lower <- 0 + CI.upper <- 1 + if (x != 0) { + CI.lower <- stats::qbeta((1 - conf.level) / 2, x, n - + x + 1) + while (acceptbin(x, n, CI.lower + tol) < (1 - + conf.level)) { + CI.lower <- CI.lower + tol + } + } + if (x != n) { + CI.upper <- stats::qbeta(1 - (1 - conf.level) / 2, x + + 1, n - x) + while (acceptbin(x, n, CI.upper - tol) < (1 - + conf.level)) { + CI.upper <- CI.upper - tol + } + } + } + ) + ci <- c(est = est, lwr.ci = max(0, CI.lower), upr.ci = min( + 1, + CI.upper + )) + if (sides == "left") { + ci[3] <- 1 + } else if (sides == "right") { + ci[2] <- 0 + } + return(ci) + } + lst <- list( + x = x, n = n, conf.level = conf.level, sides = sides, + method = method, rand = rand + ) + maxdim <- max(unlist(lapply(lst, length))) + lgp <- lapply(lst, rep, length.out = maxdim) + lgn <- h_recycle(x = if (is.null(names(x))) { + paste("x", seq_along(x), sep = ".") + } else { + names(x) + }, n = if (is.null(names(n))) { + paste("n", seq_along(n), sep = ".") + } else { + names(n) + }, conf.level = conf.level, sides = sides, method = method) + xn <- apply(as.data.frame(lgn[sapply(lgn, function(x) { + length(unique(x)) != + 1 + })]), 1, paste, collapse = ":") + res <- t(sapply(1:maxdim, function(i) { + iBinomCI( + x = lgp$x[i], + n = lgp$n[i], conf.level = lgp$conf.level[i], sides = lgp$sides[i], + method = lgp$method[i], rand = lgp$rand[i] + ) + })) + colnames(res)[1] <- c("est") + rownames(res) <- xn + return(res) +} diff --git a/R/prop_diff.R b/R/prop_diff.R index b82051b027..7c8f54236c 100644 --- a/R/prop_diff.R +++ b/R/prop_diff.R @@ -89,36 +89,44 @@ d_proportion_diff <- function(conf_level, #' set.seed(2) #' rsp <- sample(c(TRUE, FALSE), replace = TRUE, size = 20) #' grp <- factor(c(rep("A", 10), rep("B", 10))) -#' prop_diff_wald(rsp = rsp, grp = grp, conf_level = 0.90, correct = FALSE) +#' prop_diff_wald(rsp = rsp, grp = grp, conf_level = 0.95, correct = FALSE) #' #' @export prop_diff_wald <- function(rsp, grp, - conf_level, - correct) { + conf_level = 0.95, + correct = FALSE) { + if (isTRUE(correct)) { + mthd <- "waldcc" + } else { + mthd <- "wald" + } grp <- as_factor_keep_attributes(grp) check_diff_prop_ci( rsp = rsp, grp = grp, conf_level = conf_level, correct = correct ) - - p_grp <- tapply(rsp, grp, mean) diff_ci <- if (all(rsp == rsp[1])) { c(NA, NA) } else { - stats::prop.test( - table(grp, rsp), - correct = correct, - conf.level = conf_level - )$conf.int[1:2] + # check if binary response is coded as logical + checkmate::assert_logical(rsp, any.missing = FALSE) + checkmate::assert_factor(grp, len = length(rsp), any.missing = FALSE, n.levels = 2) + + tbl <- table(grp, factor(rsp, levels = c(TRUE, FALSE))) + # x1 and n1 are non-reference groups. + desctools_binom( + x1 = tbl[2], n1 = sum(tbl[2], tbl[4]), + x2 = tbl[1], n2 = sum(tbl[1], tbl[3]), + conf.level = conf_level, + method = mthd + ) } - list( - diff = unname(diff(p_grp)), - diff_ci = diff_ci + "diff" = unname(diff_ci[, "est"]), + "diff_ci" = unname(diff_ci[, c("lwr.ci", "upr.ci")]) ) } - #' @describeIn prop_diff Anderson-Hauck confidence interval. #' #' @examples @@ -140,21 +148,23 @@ prop_diff_ha <- function(rsp, grp <- as_factor_keep_attributes(grp) check_diff_prop_ci(rsp = rsp, grp = grp, conf_level = conf_level) - n_grp <- tapply(rsp, grp, length) - p_grp <- tapply(rsp, grp, mean) - diff_p <- unname(diff(p_grp)) - z <- stats::qnorm((1 + conf_level) / 2) - err <- 1 / - (2 * min(n_grp)) + z * sqrt(sum(p_grp * (1 - p_grp) / (n_grp - 1))) - l_ci <- max(-1, diff_p - err) - u_ci <- min(1, diff_p + err) + tbl <- table(grp, factor(rsp, levels = c(TRUE, FALSE))) + # x1 and n1 are non-reference groups. + ci <- desctools_binom( + x1 = tbl[2], n1 = sum(tbl[2], tbl[4]), + x2 = tbl[1], n2 = sum(tbl[1], tbl[3]), + conf.level = conf_level, + method = "ha" + ) list( - "diff" = diff_p, - "diff_ci" = c(l_ci, u_ci) + "diff" = unname(ci[, "est"]), + "diff_ci" = unname(ci[, c("lwr.ci", "upr.ci")]) ) } + + #' @describeIn prop_diff Newcombe confidence interval. It is based on #' the Wilson score confidence interval for a single binomial proportion. #' @@ -175,24 +185,27 @@ prop_diff_nc <- function(rsp, grp, conf_level, correct = FALSE) { + if (isTRUE(correct)) { + mthd <- "scorecc" + } else { + mthd <- "score" + } grp <- as_factor_keep_attributes(grp) check_diff_prop_ci(rsp = rsp, grp = grp, conf_level = conf_level) - # Source: - # https://www.lexjansen.com/wuss/2016/127_Final_Paper_PDF.pdf p_grp <- tapply(rsp, grp, mean) diff_p <- unname(diff(p_grp)) - x_grp <- split(rsp, f = grp) - ci_grp <- lapply(x_grp, FUN = prop_wilson, correct = correct, conf_level = conf_level) - l1 <- ci_grp[[1]][1] - u1 <- ci_grp[[1]][2] - l2 <- ci_grp[[2]][1] - u2 <- ci_grp[[2]][2] - l_ci <- max(-1, diff_p - sqrt((u1 - p_grp[1])^2 + (p_grp[2] - l2)^2)) - u_ci <- min(1, diff_p + sqrt((p_grp[1] - l1)^2 + (u2 - p_grp[2])^2)) + tbl <- table(grp, factor(rsp, levels = c(TRUE, FALSE))) + ci <- desctools_binom( + # x1 and n1 are non-reference groups. + x1 = tbl[2], n1 = sum(tbl[2], tbl[4]), + x2 = tbl[1], n2 = sum(tbl[1], tbl[3]), + conf.level = conf_level, + method = mthd + ) list( - "diff" = diff_p, - "diff_ci" = c(l_ci, u_ci) + "diff" = unname(ci[, "est"]), + "diff_ci" = unname(ci[, c("lwr.ci", "upr.ci")]) ) } diff --git a/R/prop_diff_test.R b/R/prop_diff_test.R index f34ce91b45..010b003ae3 100644 --- a/R/prop_diff_test.R +++ b/R/prop_diff_test.R @@ -87,10 +87,12 @@ prop_cmh <- function(ary) { #' @describeIn prop_diff_test performs the Chi-Squared test with Schouten -#' correction ([Schouten 1980]( -#' https://onlinelibrary.wiley.com/doi/abs/10.1002/bimj.4710220305)). +#' correction. #' @order 2 #' +#' @seealso For information on the Schouten correction (Schouten, 1980), +#' visit https://onlinelibrary.wiley.com/doi/abs/10.1002/bimj.4710220305. +#' #' @examples #' ## Chi-Squared test + Schouten correction. #' # Internal function - prop_schouten diff --git a/R/stat.R b/R/stat.R index d1293153ed..bdb6ccddce 100644 --- a/R/stat.R +++ b/R/stat.R @@ -136,3 +136,43 @@ stat_median_ci <- function(x, return(ci) } + +#' p-Value of the Mean +#' +#' @description `r lifecycle::badge("stable")` +#' +#' Convenient function for calculating the two-sided p-value of the mean. +#' +#' @inheritParams argument_convention +#' @param n_min (`number`)\cr a minimum number of non-missing `x` to estimate +#' the p-value of the mean. +#' @param test_mean (`number`)\cr mean value to test under the null hypothesis. +#' +#' @examples +#' stat_mean_pval(sample(10)) +#' +#' stat_mean_pval(rnorm(10), test_mean = 0.5) +#' +#' @export +stat_mean_pval <- function(x, + na.rm = TRUE, + n_min = 2, + test_mean = 0) { + if (na.rm) { + x <- stats::na.omit(x) + } + n <- length(x) + + x_mean <- mean(x) + x_sd <- stats::sd(x) + + if (n < n_min) { + pv <- c(p_value = NA_real_) + } else { + x_se <- stats::sd(x) / sqrt(n) + ttest <- (x_mean - test_mean) / x_se + pv <- c(p_value = 2 * stats::pt(-abs(ttest), df = n - 1)) + } + + return(pv) +} diff --git a/R/summarize_change.R b/R/summarize_change.R index 26b2cb9584..8aeb192a4d 100644 --- a/R/summarize_change.R +++ b/R/summarize_change.R @@ -75,13 +75,16 @@ a_change_from_baseline <- make_afun( .formats = c( n = "xx", mean_sd = "xx.xx (xx.xx)", + mean_se = "xx.xx (xx.xx)", median = "xx.xx", range = "xx.xx - xx.xx", mean_ci = "(xx.xx, xx.xx)", - median_ci = "(xx.xx, xx.xx)" + median_ci = "(xx.xx, xx.xx)", + mean_pval = "xx.xx" ), .labels = c( mean_sd = "Mean (SD)", + mean_se = "Mean (SE)", median = "Median", range = "Min - Max" ) diff --git a/R/summarize_variables.R b/R/summarize_variables.R index dddbff8eda..cd1daab909 100644 --- a/R/summarize_variables.R +++ b/R/summarize_variables.R @@ -10,18 +10,21 @@ #' @param quantile_type (`numeric`) \cr between 1 and 9 selecting quantile algorithms to be used. \cr #' Default is set to `2` as this matches the default quantile algorithm in SAS `proc univariate` set by `QNTLDEF=5`. #' This differs from R's default. See more about `type` in [stats::quantile()]. +#' @param test_mean (`numeric`) \cr to test against the mean under the null hypothesis when calculating p-value. #' #' @return A list of components with the same names as the arguments. #' @export #' control_summarize_vars <- function(conf_level = 0.95, quantiles = c(0.25, 0.75), - quantile_type = 2) { + quantile_type = 2, + test_mean = 0) { checkmate::assert_vector(quantiles, len = 2) checkmate::assert_int(quantile_type, lower = 1, upper = 9) + checkmate::assert_numeric(test_mean) nullo <- lapply(quantiles, assert_proportion_value) assert_proportion_value(conf_level) - list(conf_level = conf_level, quantiles = quantiles, quantile_type = quantile_type) + list(conf_level = conf_level, quantiles = quantiles, quantile_type = quantile_type, test_mean = test_mean) } #' Format Function for Descriptive Statistics @@ -51,9 +54,11 @@ summary_formats <- function(type = "numeric") { sd = "xx.x", se = "xx.x", mean_sd = "xx.x (xx.x)", + mean_se = "xx.x (xx.x)", mean_ci = "(xx.xx, xx.xx)", mean_sei = "(xx.xx, xx.xx)", mean_sdi = "(xx.xx, xx.xx)", + mean_pval = "xx.xx", median = "xx.x", mad = "xx.x", median_ci = "(xx.xx, xx.xx)", @@ -83,6 +88,7 @@ summary_labels <- function() { sd = "SD", se = "SE", mean_sd = "Mean (SD)", + mean_se = "Mean (SE)", median = "Median", mad = "Median Absolute Deviation", iqr = "IQR", @@ -142,6 +148,7 @@ s_summary <- function(x, #' * `quantiles`: numeric vector of length two to specify the quantiles. #' * `quantile_type` (`numeric`) \cr between 1 and 9 selecting quantile algorithms to be used. \cr #' See more about `type` in [stats::quantile()]. +#' * `test_mean`: (`numeric`) \cr to test against the mean under the null hypothesis when calculating p-value. #' #' @return If `x` is of class `numeric`, returns a list with named items: \cr #' - `n`: the [length()] of `x`. @@ -150,9 +157,11 @@ s_summary <- function(x, #' - `sd`: the [stats::sd()] of `x`. #' - `se`: the standard error of `x` mean, i.e.: (`sd()/sqrt(length())]`). #' - `mean_sd`: the [mean()] and [stats::sd()] of `x`. +#' - `mean_se`: the [mean()] of `x` and its standard error (see above). #' - `mean_ci`: the CI for the mean of `x` (from [stat_mean_ci()]). #' - `mean_sei`: the SE interval for the mean of `x`, i.e.: ([mean()] -/+ [stats::sd()]/[sqrt()]). #' - `mean_sdi`: the SD interval for the mean of `x`, i.e.: ([mean()] -/+ [stats::sd()]). +#' - `mean_pval`: the two-sided p-value of the mean of `x` (from [stat_mean_pval()]). #' - `median`: the [stats::median()] of `x`. #' - `mad`: #' the median absolute deviation of `x`, i.e.: ([stats::median()] of `xc`, where `xc` = `x` - [stats::median()]). @@ -232,6 +241,8 @@ s_summary.numeric <- function(x, # nolint y$mean_sd <- c(y$mean, "sd" = stats::sd(x, na.rm = FALSE)) + y$mean_se <- c(y$mean, y$se) + mean_ci <- stat_mean_ci(x, conf_level = control$conf_level, na.rm = FALSE, gg_helper = FALSE) y$mean_ci <- formatters::with_label(mean_ci, paste("Mean", f_conf_level(control$conf_level))) @@ -243,6 +254,9 @@ s_summary.numeric <- function(x, # nolint names(mean_sdi) <- c("mean_sdi_lwr", "mean_sdi_upr") y$mean_sdi <- formatters::with_label(mean_sdi, "Mean -/+ 1xSD") + mean_pval <- stat_mean_pval(x, test_mean = control$test_mean, na.rm = FALSE, n_min = 2) + y$mean_pval <- formatters::with_label(mean_pval, paste("Mean", f_pval(control$test_mean))) + y$median <- c("median" = stats::median(x, na.rm = FALSE)) y$mad <- c("mad" = stats::median(x - y$median, na.rm = FALSE)) diff --git a/R/survival_coxph_pairwise.R b/R/survival_coxph_pairwise.R index 6de8e63b86..592c5e6101 100644 --- a/R/survival_coxph_pairwise.R +++ b/R/survival_coxph_pairwise.R @@ -16,6 +16,7 @@ #' can also be set to "breslow" or "exact". See more in [survival::coxph()] #' * `conf_level`: (`proportion`)\cr confidence level of the interval for HR. #' +#' @importFrom stats pchisq #' @name survival_coxph_pairwise #' NULL diff --git a/R/utils.R b/R/utils.R index 50aa0e839c..c190e4d353 100644 --- a/R/utils.R +++ b/R/utils.R @@ -55,6 +55,19 @@ f_conf_level <- function(conf_level) { paste0(conf_level * 100, "% CI") } +#' Utility function to create label for p-value +#' +#' @description `r lifecycle::badge("stable")` +#' +#' @param test_mean (`number`)\cr mean value to test under the null hypothesis. +#' @return a `string` +#' +#' @export +f_pval <- function(test_mean) { + checkmate::assert_numeric(test_mean, len = 1) + paste0("p-value (H0: mean = ", test_mean, ")") +} + #' Utility function to return a named list of covariate names. #' #' @param covariates (`character`)\cr a vector that can contain single variable names (such as @@ -222,12 +235,9 @@ empty_vector_if_na <- function(x) { #' @return A `list` where each element combines corresponding elements of `x` and `y`. #' #' @examples -#' # Internal function - combine_vectors -#' \dontrun{ #' combine_vectors(1:3, 4:6) -#' } #' -#' @keywords internal +#' @export combine_vectors <- function(x, y) { checkmate::assert_vector(x) checkmate::assert_vector(y, len = length(x)) diff --git a/_pkgdown.yml b/_pkgdown.yml index 3c7de91d47..eef6bd1cec 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -1,3 +1,4 @@ +--- url: https://insightsengineering.github.io/tern template: @@ -16,8 +17,9 @@ reference: - kaplan_meier - title: Control Functions - desc: Control functions capture options in lists, and take care of defaults (and if it makes sense also checks). - They avoid cluttering of function signatures with long lists of single arguments. + desc: Control functions capture options in lists, and take care of defaults + (and if it makes sense also checks). They avoid cluttering of + function signatures with long lists of single arguments. contents: - control_coxph - control_coxreg @@ -30,9 +32,10 @@ reference: - control_surv_timepoint - title: Statistics Functions - desc: Statistics functions should do the computation of the numbers that are tabulated later. - In order to separate computation from formatting, they should not take care of `rcell` type - formatting themselves. + desc: Statistics functions should do the computation of the numbers that + are tabulated later. In order to separate computation from + formatting, they should not take care of `rcell` type formatting + themselves. contents: - s_compare - s_count_cumulative @@ -55,59 +58,60 @@ reference: - s_summary - title: Formatted Analysis functions - desc: These have the same arguments as the corresponding statistics functions, and can be further - customized by calling `rtables::make_afun()` on them. They are used as `afun` in + desc: These have the same arguments as the corresponding statistics + functions, and can be further customized by calling + `rtables::make_afun()` on them. They are used as `afun` in `rtables::analyze()`. contents: - - a_count_cumulative - - a_count_missed_doses - - a_count_occurrences - - a_count_occurrences_by_grade - - a_count_patients_with_event - - a_count_patients_with_flags - - a_odds_ratio - - a_proportion - - create_afun_compare - - create_afun_summary + - a_count_cumulative + - a_count_missed_doses + - a_count_occurrences + - a_count_occurrences_by_grade + - a_count_patients_with_event + - a_count_patients_with_flags + - a_odds_ratio + - a_proportion + - create_afun_compare + - create_afun_summary - title: Analyze Functions - desc: Analyze Functions are used in combination with the rtables layout functions, in the pipeline which creates - the table. + desc: Analyze Functions are used in combination with the rtables layout + functions, in the pipeline which creates the table. contents: - - compare_vars - - count_abnormal - - count_abnormal_by_baseline - - count_abnormal_by_marked - - count_abnormal_by_worst_grade - - count_abnormal_lab_worsen_by_baseline - - count_cumulative - - count_missed_doses - - count_occurrences - - count_occurrences_by_grade - - count_patients_with_event - - count_patients_with_flags - - count_values_funs - - estimate_incidence_rate - - estimate_multinomial_response - - estimate_odds_ratio - - estimate_proportion - - estimate_proportion_diff - - response_biomarkers_subgroups - - summarize_ancova - - summarize_change - - summarize_colvars - - summarize_coxreg - - summarize_logistic - - summarize_num_patients - - summarize_occurrences_by_grade - - summarize_patients_events_in_cols - - summarize_patients_exposure_in_cols - - summarize_vars_in_cols - - summarize_vars - - tabulate_rsp_subgroups - - tabulate_survival_biomarkers - - tabulate_survival_subgroups - - test_proportion_diff + - compare_vars + - count_abnormal + - count_abnormal_by_baseline + - count_abnormal_by_marked + - count_abnormal_by_worst_grade + - count_abnormal_lab_worsen_by_baseline + - count_cumulative + - count_missed_doses + - count_occurrences + - count_occurrences_by_grade + - count_patients_with_event + - count_patients_with_flags + - count_values_funs + - estimate_incidence_rate + - estimate_multinomial_response + - estimate_odds_ratio + - estimate_proportion + - estimate_proportion_diff + - response_biomarkers_subgroups + - summarize_ancova + - summarize_change + - summarize_colvars + - summarize_coxreg + - summarize_logistic + - summarize_num_patients + - summarize_occurrences_by_grade + - summarize_patients_events_in_cols + - summarize_patients_exposure_in_cols + - summarize_vars_in_cols + - summarize_vars + - tabulate_rsp_subgroups + - tabulate_survival_biomarkers + - tabulate_survival_subgroups + - test_proportion_diff - title: Analysis Helper Functions desc: These functions are useful to help definining the analysis. @@ -162,6 +166,7 @@ reference: - h_pkparam_sort - h_proportion_df - h_proportion_subgroups_df + - h_recycle - h_response_biomarkers_subgroups - h_rsp_to_logistic_variables - h_simple_term_labels @@ -189,9 +194,12 @@ reference: - prop_diff_ha - prop_diff_nc - prop_diff_wald + - desctools_binom + - desctools_binomci - title: Model Specific Functions - desc: These functions help with fitting or extracting results from specific models. + desc: These functions help with fitting or extracting results from specific + models. contents: - estimate_coef - extract_rsp_subgroups @@ -224,13 +232,15 @@ reference: - g_waterfall - title: rtables Helper Functions - desc: These functions help to work with the new rtables package and might be moved there later. + desc: These functions help to work with the new rtables package and might + be moved there later. contents: - add_rowcounts - append_varlabels - as.rtable - as.rtable.data.frame - combine_groups + - combine_vectors - h_col_counts - h_col_indices - h_content_first_row @@ -241,7 +251,8 @@ reference: - to_string_matrix - title: rtables Formatting Functions - desc: These functions provide customized formatting rules to work with the rtables package. + desc: These functions provide customized formatting rules to work with the + rtables package. contents: - format_count_fraction - format_extreme_values @@ -285,6 +296,7 @@ reference: - h_xticks - stack_grobs - stat_mean_ci + - stat_mean_pval - stat_median_ci - title: Data Helper Functions @@ -313,6 +325,7 @@ reference: - fct_discard - fct_explicit_na_if - f_conf_level + - f_pval - h_data_plot - sas_na - to_n diff --git a/inst/WORDLIST b/inst/WORDLIST index 87828b2ecc..278bac0238 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -24,23 +24,36 @@ Pharmacokinetics Sabanés Satterthwaite Schouten -timepoint TLG TLGs Wojciak Wojciech Yuyao +agresti +arcsine +bimj biomarker biomarkers +blaker breslow +clopper coercible +coull coxph +doi efron funder +https +jeffreys +lik lineplot -mmrm +matric +midp nd +onlinelibrary paramcd +pearson +pratt pre responder responders @@ -49,8 +62,12 @@ songy subtable subtables surv +timepoint unstratified waddella wald +wiley +wilson +wilsoncc wojciak wojciech diff --git a/man/combine_vectors.Rd b/man/combine_vectors.Rd index 528d530fef..4b5ff87a82 100644 --- a/man/combine_vectors.Rd +++ b/man/combine_vectors.Rd @@ -18,10 +18,6 @@ A \code{list} where each element combines corresponding elements of \code{x} and Combine Two Vectors Element Wise } \examples{ -# Internal function - combine_vectors -\dontrun{ combine_vectors(1:3, 4:6) -} } -\keyword{internal} diff --git a/man/control_summarize_vars.Rd b/man/control_summarize_vars.Rd index 2984d0a700..230839741a 100644 --- a/man/control_summarize_vars.Rd +++ b/man/control_summarize_vars.Rd @@ -7,7 +7,8 @@ control_summarize_vars( conf_level = 0.95, quantiles = c(0.25, 0.75), - quantile_type = 2 + quantile_type = 2, + test_mean = 0 ) } \arguments{ @@ -18,6 +19,8 @@ control_summarize_vars( \item{quantile_type}{(\code{numeric}) \cr between 1 and 9 selecting quantile algorithms to be used. \cr Default is set to \code{2} as this matches the default quantile algorithm in SAS \verb{proc univariate} set by \code{QNTLDEF=5}. This differs from R's default. See more about \code{type} in \code{\link[stats:quantile]{stats::quantile()}}.} + +\item{test_mean}{(\code{numeric}) \cr to test against the mean under the null hypothesis when calculating p-value.} } \value{ A list of components with the same names as the arguments. diff --git a/man/desctools_binom.Rd b/man/desctools_binom.Rd new file mode 100644 index 0000000000..c9663bd4db --- /dev/null +++ b/man/desctools_binom.Rd @@ -0,0 +1,91 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/desctools_binom_diff.R +\name{desctools_binom} +\alias{desctools_binom} +\alias{h_recycle} +\alias{desctools_binomci} +\title{Confidence Intervals for a Difference of Binomials} +\usage{ +h_recycle(...) + +desctools_binom( + x1, + n1, + x2, + n2, + conf.level = 0.95, + sides = c("two.sided", "left", "right"), + method = c("ac", "wald", "waldcc", "score", "scorecc", "mn", "mee", "blj", "ha", "hal", + "jp") +) + +desctools_binomci( + x, + n, + conf.level = 0.95, + sides = c("two.sided", "left", "right"), + method = c("wilson", "wald", "waldcc", "agresti-coull", "jeffreys", "modified wilson", + "wilsoncc", "modified jeffreys", "clopper-pearson", "arcsine", "logit", "witting", + "pratt", "midp", "lik", "blaker"), + rand = 123, + tol = 1e-05 +) +} +\arguments{ +\item{conf.level}{\cr confidence level, defaults to 0.95} + +\item{sides}{\cr a character string specifying the side of the confidence interval. Must be one of "two-sided" (default), "left" or "right".} + +\item{method}{\cr character string specifying which method to use. Can be one out of: "wald", "wilson", "wilsoncc", "agresti-coull", "jeffreys", "modified wilson", "modified jeffreys", "clopper-pearson", "arcsine. +", "logit", "witting", "pratt", "midp", "lik" and "blaker" + +@return A matric with 3 columns containing: +\itemize{ +\item \code{est}: estimate of proportion difference. +\item \code{lwrci}: lower end of the confidence interval +\item \code{upci}: upper end of the confidence interval. +}} + +\item{x}{\cr number of successes} + +\item{n}{\cr number of trials} + +\item{grp}{(\code{factor})\cr +vector assigning observations to one out of two groups +(e.g. reference and treatment group).} +} +\value{ +A named list of 3 values: +\itemize{ +\item \code{est}: estimate of proportion difference. +\item \code{lwrci}: estimate of lower end of the confidence interval +\item \code{upci}: estimate of upper end of the confidence interval. +} +} +\description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} + +Several confidence intervals for the difference between proportions. +} +\section{Functions}{ +\itemize{ +\item \code{h_recycle()}: This function recycles all supplied elements to the maximal dimension. + +\item \code{desctools_binom()}: Several Confidence Intervals for the difference between proportions. + +\item \code{desctools_binomci()}: Compute confidence intervals for binomial proportions. + +}} +\examples{ +# Internal function - +\dontrun{ + +set.seed(2) +rsp <- sample(c(TRUE, FALSE), replace = TRUE, size = 20) +grp <- factor(c(rep("A", 10), rep("B", 10))) +tbl <- table(grp, factor(rsp, levels = c(TRUE, FALSE))) +desctools_binom(tbl[1], sum(tbl[1], tbl[3]), tbl[2], sum(tbl[2], tbl[4]), conf.level = 0.90, method = "waldcc") +} + +} +\keyword{internal} diff --git a/man/f_pval.Rd b/man/f_pval.Rd new file mode 100644 index 0000000000..1c5f82bcb4 --- /dev/null +++ b/man/f_pval.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils.R +\name{f_pval} +\alias{f_pval} +\title{Utility function to create label for p-value} +\usage{ +f_pval(test_mean) +} +\arguments{ +\item{test_mean}{(\code{number})\cr mean value to test under the null hypothesis.} +} +\value{ +a \code{string} +} +\description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#stable}{\figure{lifecycle-stable.svg}{options: alt='[Stable]'}}}{\strong{[Stable]}} +} diff --git a/man/prop_diff.Rd b/man/prop_diff.Rd index 12694e8c62..3d8e0aaf68 100644 --- a/man/prop_diff.Rd +++ b/man/prop_diff.Rd @@ -14,7 +14,7 @@ \usage{ d_proportion_diff(conf_level, method, long = FALSE) -prop_diff_wald(rsp, grp, conf_level, correct) +prop_diff_wald(rsp, grp, conf_level = 0.95, correct = FALSE) prop_diff_ha(rsp, grp, conf_level) @@ -153,7 +153,7 @@ should be accounted for. set.seed(2) rsp <- sample(c(TRUE, FALSE), replace = TRUE, size = 20) grp <- factor(c(rep("A", 10), rep("B", 10))) -prop_diff_wald(rsp = rsp, grp = grp, conf_level = 0.90, correct = FALSE) +prop_diff_wald(rsp = rsp, grp = grp, conf_level = 0.95, correct = FALSE) # Anderson-Hauck confidence interval ## "Mid" case: 3/4 respond in group A, 1/2 respond in group B. diff --git a/man/prop_diff_test.Rd b/man/prop_diff_test.Rd index 82d09f426b..6f90991e9c 100644 --- a/man/prop_diff_test.Rd +++ b/man/prop_diff_test.Rd @@ -1,220 +1,224 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/prop_diff_test.R -\name{prop_diff_test} -\alias{prop_diff_test} -\alias{prop_chisq} -\alias{prop_schouten} -\alias{prop_fisher} -\alias{prop_cmh} -\alias{s_test_proportion_diff} -\alias{a_test_proportion_diff} -\alias{test_proportion_diff} -\title{Difference Test for Two Proportions} -\usage{ -prop_chisq(tbl) - -prop_schouten(tbl) - -prop_fisher(tbl) - -prop_cmh(ary) - -s_test_proportion_diff( - df, - .var, - .ref_group, - .in_ref_col, - variables = list(strata = NULL), - method = c("chisq", "schouten", "fisher", "cmh") -) - -a_test_proportion_diff( - df, - .var, - .ref_group, - .in_ref_col, - variables = list(strata = NULL), - method = c("chisq", "schouten", "fisher", "cmh") -) - -test_proportion_diff( - lyt, - vars, - ..., - show_labels = "hidden", - table_names = vars, - .stats = NULL, - .formats = NULL, - .labels = NULL, - .indent_mods = NULL -) -} -\arguments{ -\item{tbl}{(\code{matrix})\cr -with two groups in rows and the binary response (\code{TRUE}/\code{FALSE}) in -columns.} - -\item{ary}{(\code{array}, 3 dimensions)\cr -with two groups in rows, the binary response (\code{TRUE}/\code{FALSE}) in -columns, the strata in the third dimension.} - -\item{df}{(\verb{data frame})\cr data set containing all analysis variables.} - -\item{.var}{(\code{string})\cr single variable name that is passed by \code{rtables} when requested -by a statistics function.} - -\item{.ref_group}{(\verb{data frame} or \code{vector})\cr the data corresponding to the reference group.} - -\item{.in_ref_col}{(\code{logical})\cr \code{TRUE} when working with the reference level, \code{FALSE} otherwise.} - -\item{variables}{(named \code{list} of \code{string})\cr list of additional analysis variables.} - -\item{method}{(\code{string})\cr -one of (\code{chisq}, \code{cmh}, \code{fisher}, \code{schouten}; specifies the test used -to calculate the p-value.} - -\item{lyt}{(\code{layout})\cr input layout where analyses will be added to.} - -\item{vars}{(\code{character})\cr variable names for the primary analysis variable to be iterated over.} - -\item{...}{other arguments are passed to \code{\link[=s_test_proportion_diff]{s_test_proportion_diff()}}.} - -\item{show_labels}{label visibility: one of "default", "visible" and "hidden".} - -\item{table_names}{(\code{character})\cr this can be customized in case that the same \code{vars} are analyzed multiple times, -to avoid warnings from \code{rtables}.} - -\item{.stats}{(\code{character})\cr statistics to select for the table.} - -\item{.formats}{(named \code{character} or \code{list})\cr formats for the statistics.} - -\item{.labels}{(named \code{character})\cr labels for the statistics (without indent).} - -\item{.indent_mods}{(named \code{integer})\cr indent modifiers for the labels.} -} -\value{ -Named \code{list} with a single item \code{pval} with an attribute \code{label} -describing the method used. The p-value tests the null hypothesis that -proportions in two groups are the same. -} -\description{ -\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#stable}{\figure{lifecycle-stable.svg}{options: alt='[Stable]'}}}{\strong{[Stable]}} - -Various tests were implemented to test the difference between two -proportions. -} -\section{Functions}{ -\itemize{ -\item \code{prop_chisq()}: performs Chi-Squared test. -Internally calls \code{\link[stats:prop.test]{stats::prop.test()}}. - -\item \code{prop_schouten()}: performs the Chi-Squared test with Schouten -correction (\href{https://onlinelibrary.wiley.com/doi/abs/10.1002/bimj.4710220305}{Schouten 1980}). - -\item \code{prop_fisher()}: performs the Fisher's exact test. -Internally calls \code{\link[stats:fisher.test]{stats::fisher.test()}}. - -\item \code{prop_cmh()}: performs stratified Cochran-Mantel-Haenszel test. -Internally calls \code{\link[stats:mantelhaen.test]{stats::mantelhaen.test()}}. Note that strata with less than two observations -are automatically discarded. - -\item \code{s_test_proportion_diff()}: Statistics function which test the difference -between two proportions. - -\item \code{a_test_proportion_diff()}: Formatted Analysis function which can be further customized by calling -\code{\link[rtables:make_afun]{rtables::make_afun()}} on it. It is used as \code{afun} in \code{\link[rtables:analyze]{rtables::analyze()}}. - -\item \code{test_proportion_diff()}: Layout creating function which can be used for -creating tables, which can take statistics function arguments and -additional format arguments. - -}} -\examples{ -# Non-stratified proportion difference test - -## Data -A <- 20 -B <- 20 -set.seed(1) -rsp <- c( - sample(c(TRUE, FALSE), size = A, prob = c(3 / 4, 1 / 4), replace = TRUE), - sample(c(TRUE, FALSE), size = A, prob = c(1 / 2, 1 / 2), replace = TRUE) -) -grp <- c(rep("A", A), rep("B", B)) -tbl <- table(grp, rsp) - -## Chi-Squared test -# Internal function - prop_chisq -\dontrun{ -prop_chisq(tbl) -} - -## Chi-Squared test + Schouten correction. -# Internal function - prop_schouten -\dontrun{ -prop_schouten(tbl) -} - -## Fisher's exact test -# Internal function - prop_fisher -\dontrun{ -prop_fisher(tbl) -} - -# Stratified proportion difference test - -## Data -rsp <- sample(c(TRUE, FALSE), 100, TRUE) -grp <- factor(rep(c("A", "B"), each = 50)) -strata <- factor(rep(c("V", "W", "X", "Y", "Z"), each = 20)) -tbl <- table(grp, rsp, strata) - -## Cochran-Mantel-Haenszel test -# Internal function - prop_cmh -\dontrun{ -prop_cmh(tbl) -} - - -# Statistics function -dta <- data.frame( - rsp = sample(c(TRUE, FALSE), 100, TRUE), - grp = factor(rep(c("A", "B"), each = 50)), - strat = factor(rep(c("V", "W", "X", "Y", "Z"), each = 20)) -) - -# Internal function - s_test_proportion_diff -\dontrun{ -s_test_proportion_diff( - df = subset(dta, grp == "A"), - .var = "rsp", - .ref_group = subset(dta, grp == "B"), - .in_ref_col = FALSE, - variables = list(strata = "strat"), - method = "cmh" -) -} - -# Internal function - a_test_proportion_diff -\dontrun{ -a_test_proportion_diff( - df = subset(dta, grp == "A"), - .var = "rsp", - .ref_group = subset(dta, grp == "B"), - .in_ref_col = FALSE, - variables = list(strata = "strat"), - method = "cmh" -) -} - -# With `rtables` pipelines. -l <- basic_table() \%>\% - split_cols_by(var = "grp", ref_group = "B") \%>\% - test_proportion_diff( - vars = "rsp", - method = "cmh", variables = list(strata = "strat") - ) - -build_table(l, df = dta) -} -\keyword{internal} +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/prop_diff_test.R +\name{prop_diff_test} +\alias{prop_diff_test} +\alias{prop_chisq} +\alias{prop_schouten} +\alias{prop_fisher} +\alias{prop_cmh} +\alias{s_test_proportion_diff} +\alias{a_test_proportion_diff} +\alias{test_proportion_diff} +\title{Difference Test for Two Proportions} +\usage{ +prop_chisq(tbl) + +prop_schouten(tbl) + +prop_fisher(tbl) + +prop_cmh(ary) + +s_test_proportion_diff( + df, + .var, + .ref_group, + .in_ref_col, + variables = list(strata = NULL), + method = c("chisq", "schouten", "fisher", "cmh") +) + +a_test_proportion_diff( + df, + .var, + .ref_group, + .in_ref_col, + variables = list(strata = NULL), + method = c("chisq", "schouten", "fisher", "cmh") +) + +test_proportion_diff( + lyt, + vars, + ..., + show_labels = "hidden", + table_names = vars, + .stats = NULL, + .formats = NULL, + .labels = NULL, + .indent_mods = NULL +) +} +\arguments{ +\item{tbl}{(\code{matrix})\cr +with two groups in rows and the binary response (\code{TRUE}/\code{FALSE}) in +columns.} + +\item{ary}{(\code{array}, 3 dimensions)\cr +with two groups in rows, the binary response (\code{TRUE}/\code{FALSE}) in +columns, the strata in the third dimension.} + +\item{df}{(\verb{data frame})\cr data set containing all analysis variables.} + +\item{.var}{(\code{string})\cr single variable name that is passed by \code{rtables} when requested +by a statistics function.} + +\item{.ref_group}{(\verb{data frame} or \code{vector})\cr the data corresponding to the reference group.} + +\item{.in_ref_col}{(\code{logical})\cr \code{TRUE} when working with the reference level, \code{FALSE} otherwise.} + +\item{variables}{(named \code{list} of \code{string})\cr list of additional analysis variables.} + +\item{method}{(\code{string})\cr +one of (\code{chisq}, \code{cmh}, \code{fisher}, \code{schouten}; specifies the test used +to calculate the p-value.} + +\item{lyt}{(\code{layout})\cr input layout where analyses will be added to.} + +\item{vars}{(\code{character})\cr variable names for the primary analysis variable to be iterated over.} + +\item{...}{other arguments are passed to \code{\link[=s_test_proportion_diff]{s_test_proportion_diff()}}.} + +\item{show_labels}{label visibility: one of "default", "visible" and "hidden".} + +\item{table_names}{(\code{character})\cr this can be customized in case that the same \code{vars} are analyzed multiple times, +to avoid warnings from \code{rtables}.} + +\item{.stats}{(\code{character})\cr statistics to select for the table.} + +\item{.formats}{(named \code{character} or \code{list})\cr formats for the statistics.} + +\item{.labels}{(named \code{character})\cr labels for the statistics (without indent).} + +\item{.indent_mods}{(named \code{integer})\cr indent modifiers for the labels.} +} +\value{ +Named \code{list} with a single item \code{pval} with an attribute \code{label} +describing the method used. The p-value tests the null hypothesis that +proportions in two groups are the same. +} +\description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#stable}{\figure{lifecycle-stable.svg}{options: alt='[Stable]'}}}{\strong{[Stable]}} + +Various tests were implemented to test the difference between two +proportions. +} +\section{Functions}{ +\itemize{ +\item \code{prop_chisq()}: performs Chi-Squared test. +Internally calls \code{\link[stats:prop.test]{stats::prop.test()}}. + +\item \code{prop_schouten()}: performs the Chi-Squared test with Schouten +correction. + +\item \code{prop_fisher()}: performs the Fisher's exact test. +Internally calls \code{\link[stats:fisher.test]{stats::fisher.test()}}. + +\item \code{prop_cmh()}: performs stratified Cochran-Mantel-Haenszel test. +Internally calls \code{\link[stats:mantelhaen.test]{stats::mantelhaen.test()}}. Note that strata with less than two observations +are automatically discarded. + +\item \code{s_test_proportion_diff()}: Statistics function which test the difference +between two proportions. + +\item \code{a_test_proportion_diff()}: Formatted Analysis function which can be further customized by calling +\code{\link[rtables:make_afun]{rtables::make_afun()}} on it. It is used as \code{afun} in \code{\link[rtables:analyze]{rtables::analyze()}}. + +\item \code{test_proportion_diff()}: Layout creating function which can be used for +creating tables, which can take statistics function arguments and +additional format arguments. + +}} +\examples{ +# Non-stratified proportion difference test + +## Data +A <- 20 +B <- 20 +set.seed(1) +rsp <- c( + sample(c(TRUE, FALSE), size = A, prob = c(3 / 4, 1 / 4), replace = TRUE), + sample(c(TRUE, FALSE), size = A, prob = c(1 / 2, 1 / 2), replace = TRUE) +) +grp <- c(rep("A", A), rep("B", B)) +tbl <- table(grp, rsp) + +## Chi-Squared test +# Internal function - prop_chisq +\dontrun{ +prop_chisq(tbl) +} + +## Chi-Squared test + Schouten correction. +# Internal function - prop_schouten +\dontrun{ +prop_schouten(tbl) +} + +## Fisher's exact test +# Internal function - prop_fisher +\dontrun{ +prop_fisher(tbl) +} + +# Stratified proportion difference test + +## Data +rsp <- sample(c(TRUE, FALSE), 100, TRUE) +grp <- factor(rep(c("A", "B"), each = 50)) +strata <- factor(rep(c("V", "W", "X", "Y", "Z"), each = 20)) +tbl <- table(grp, rsp, strata) + +## Cochran-Mantel-Haenszel test +# Internal function - prop_cmh +\dontrun{ +prop_cmh(tbl) +} + + +# Statistics function +dta <- data.frame( + rsp = sample(c(TRUE, FALSE), 100, TRUE), + grp = factor(rep(c("A", "B"), each = 50)), + strat = factor(rep(c("V", "W", "X", "Y", "Z"), each = 20)) +) + +# Internal function - s_test_proportion_diff +\dontrun{ +s_test_proportion_diff( + df = subset(dta, grp == "A"), + .var = "rsp", + .ref_group = subset(dta, grp == "B"), + .in_ref_col = FALSE, + variables = list(strata = "strat"), + method = "cmh" +) +} + +# Internal function - a_test_proportion_diff +\dontrun{ +a_test_proportion_diff( + df = subset(dta, grp == "A"), + .var = "rsp", + .ref_group = subset(dta, grp == "B"), + .in_ref_col = FALSE, + variables = list(strata = "strat"), + method = "cmh" +) +} + +# With `rtables` pipelines. +l <- basic_table() \%>\% + split_cols_by(var = "grp", ref_group = "B") \%>\% + test_proportion_diff( + vars = "rsp", + method = "cmh", variables = list(strata = "strat") + ) + +build_table(l, df = dta) +} +\seealso{ +For information on the Schouten correction (Schouten, 1980), +visit https://onlinelibrary.wiley.com/doi/abs/10.1002/bimj.4710220305. +} +\keyword{internal} diff --git a/man/stat_mean_pval.Rd b/man/stat_mean_pval.Rd new file mode 100644 index 0000000000..49f7ae4380 --- /dev/null +++ b/man/stat_mean_pval.Rd @@ -0,0 +1,29 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stat.R +\name{stat_mean_pval} +\alias{stat_mean_pval} +\title{p-Value of the Mean} +\usage{ +stat_mean_pval(x, na.rm = TRUE, n_min = 2, test_mean = 0) +} +\arguments{ +\item{x}{(\code{numeric})\cr vector of numbers we want to analyze.} + +\item{na.rm}{(\code{flag})\cr whether \code{NA} values should be removed from \code{x} prior to analysis.} + +\item{n_min}{(\code{number})\cr a minimum number of non-missing \code{x} to estimate +the p-value of the mean.} + +\item{test_mean}{(\code{number})\cr mean value to test under the null hypothesis.} +} +\description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#stable}{\figure{lifecycle-stable.svg}{options: alt='[Stable]'}}}{\strong{[Stable]}} + +Convenient function for calculating the two-sided p-value of the mean. +} +\examples{ +stat_mean_pval(sample(10)) + +stat_mean_pval(rnorm(10), test_mean = 0.5) + +} diff --git a/man/summarize_variables.Rd b/man/summarize_variables.Rd index 188f59710d..958d53c698 100644 --- a/man/summarize_variables.Rd +++ b/man/summarize_variables.Rd @@ -151,6 +151,7 @@ the helper function \code{\link[=control_summarize_vars]{control_summarize_vars( \item \code{quantiles}: numeric vector of length two to specify the quantiles. \item \code{quantile_type} (\code{numeric}) \cr between 1 and 9 selecting quantile algorithms to be used. \cr See more about \code{type} in \code{\link[stats:quantile]{stats::quantile()}}. +\item \code{test_mean}: (\code{numeric}) \cr to test against the mean under the null hypothesis when calculating p-value. }} \item{verbose}{defaults to \code{TRUE}. It prints out warnings and messages. It is mainly used @@ -186,9 +187,11 @@ If \code{x} is of class \code{numeric}, returns a list with named items: \cr \item \code{sd}: the \code{\link[stats:sd]{stats::sd()}} of \code{x}. \item \code{se}: the standard error of \code{x} mean, i.e.: (\verb{sd()/sqrt(length())]}). \item \code{mean_sd}: the \code{\link[=mean]{mean()}} and \code{\link[stats:sd]{stats::sd()}} of \code{x}. +\item \code{mean_se}: the \code{\link[=mean]{mean()}} of \code{x} and its standard error (see above). \item \code{mean_ci}: the CI for the mean of \code{x} (from \code{\link[=stat_mean_ci]{stat_mean_ci()}}). \item \code{mean_sei}: the SE interval for the mean of \code{x}, i.e.: (\code{\link[=mean]{mean()}} -/+ \code{\link[stats:sd]{stats::sd()}}/\code{\link[=sqrt]{sqrt()}}). \item \code{mean_sdi}: the SD interval for the mean of \code{x}, i.e.: (\code{\link[=mean]{mean()}} -/+ \code{\link[stats:sd]{stats::sd()}}). +\item \code{mean_pval}: the two-sided p-value of the mean of \code{x} (from \code{\link[=stat_mean_pval]{stat_mean_pval()}}). \item \code{median}: the \code{\link[stats:median]{stats::median()}} of \code{x}. \item \code{mad}: the median absolute deviation of \code{x}, i.e.: (\code{\link[stats:median]{stats::median()}} of \code{xc}, where \code{xc} = \code{x} - \code{\link[stats:median]{stats::median()}}). diff --git a/tests/testthat/test-compare_variables.R b/tests/testthat/test-compare_variables.R index 19ade7bf99..910bcb5adf 100644 --- a/tests/testthat/test-compare_variables.R +++ b/tests/testthat/test-compare_variables.R @@ -5,9 +5,9 @@ testthat::test_that("s_compare works for numeric", { .in_ref_col = FALSE )) testthat::expect_named(result, c( - "n", "sum", "mean", "sd", "se", "mean_sd", "mean_ci", "mean_sei", "mean_sdi", - "median", "mad", "median_ci", "quantiles", "iqr", "range", - "min", "max", "cv", "geom_mean", "geom_mean_ci", "geom_cv", "pval" + "n", "sum", "mean", "sd", "se", "mean_sd", "mean_se", "mean_ci", "mean_sei", + "mean_sdi", "mean_pval", "median", "mad", "median_ci", "quantiles", "iqr", + "range", "min", "max", "cv", "geom_mean", "geom_mean_ci", "geom_cv", "pval" )) }) diff --git a/tests/testthat/test-count_patients_with_event.R b/tests/testthat/test-count_patients_with_event.R index 3e3282f082..f006d12c9e 100644 --- a/tests/testthat/test-count_patients_with_event.R +++ b/tests/testthat/test-count_patients_with_event.R @@ -237,8 +237,8 @@ testthat::test_that("count_patients_with_flags works as expected when specifying }) testthat::test_that("count_patients_with_flags works with label row specified", { - adsl <- scda::synthetic_cdisc_data("rcd_2021_07_07")$adsl - adae <- scda::synthetic_cdisc_data("rcd_2021_07_07")$adae + adsl <- scda::synthetic_cdisc_data("rcd_2022_02_28")$adsl + adae <- scda::synthetic_cdisc_data("rcd_2022_02_28")$adae # Create custom flags: adae <- adae %>% diff --git a/tests/testthat/test-h_response_biomarkers_subgroups.R b/tests/testthat/test-h_response_biomarkers_subgroups.R index dac16ae11c..b6d477a6a5 100644 --- a/tests/testthat/test-h_response_biomarkers_subgroups.R +++ b/tests/testthat/test-h_response_biomarkers_subgroups.R @@ -17,7 +17,7 @@ preprocess_adrs <- function(adrs) { ) } -adrs <- synthetic_cdisc_data("rcd_2021_05_05")$adrs +adrs <- synthetic_cdisc_data("rcd_2022_02_28")$adrs # h_rsp_to_logistic_variables ---- @@ -56,7 +56,7 @@ testthat::test_that("h_logistic_mult_cont_df works as expected", { )) expected <- data.frame( biomarker = c("BMRKR1", "AGE"), - biomarker_label = c("Continous Level Biomarker 1", "Age"), + biomarker_label = c("Continuous Level Biomarker 1", "Age"), n_tot = c(400L, 400L), n_rsp = c(336L, 336L), prop = c(0.84, 0.84), @@ -85,7 +85,7 @@ testthat::test_that("h_logistic_mult_cont_df returns missing values if data is e )) expected <- data.frame( biomarker = c("BMRKR1", "AGE"), - biomarker_label = c("Continous Level Biomarker 1", "Age"), + biomarker_label = c("Continuous Level Biomarker 1", "Age"), n_tot = c(0L, 0L), n_rsp = c(0L, 0L), prop = c(NA, NA), @@ -115,7 +115,7 @@ testthat::test_that("h_logistic_mult_cont_df also works with response not being )) expected <- data.frame( biomarker = c("BMRKR1", "AGE"), - biomarker_label = c("Continous Level Biomarker 1", "Age"), + biomarker_label = c("Continuous Level Biomarker 1", "Age"), n_tot = c(400L, 400L), n_rsp = c(336L, 336L), prop = c(0.84, 0.84), diff --git a/tests/testthat/test-h_survival_biomarkers_subgroups.R b/tests/testthat/test-h_survival_biomarkers_subgroups.R index ee78a711ca..416de853e7 100644 --- a/tests/testthat/test-h_survival_biomarkers_subgroups.R +++ b/tests/testthat/test-h_survival_biomarkers_subgroups.R @@ -20,7 +20,7 @@ preprocess_adtte <- function(adtte) { ) } -adtte <- synthetic_cdisc_data("rcd_2021_05_05")$adtte +adtte <- synthetic_cdisc_data("rcd_2022_02_28")$adtte # h_surv_to_coxreg_variables ---- @@ -62,7 +62,7 @@ testthat::test_that("h_coxreg_mult_cont_df works as expected", { )) expected <- data.frame( biomarker = c("BMRKR1", "AGE"), - biomarker_label = c("Continous Level Biomarker 1", "Age"), + biomarker_label = c("Continuous Level Biomarker 1", "Age"), n_tot = c(400L, 400L), n_tot_events = c(282, 282), median = c(680.959764183532, 680.959764183532), @@ -92,7 +92,7 @@ testthat::test_that("h_coxreg_mult_cont_df returns missing values if data is emp )) expected <- data.frame( biomarker = c("BMRKR1", "AGE"), - biomarker_label = c("Continous Level Biomarker 1", "Age"), + biomarker_label = c("Continuous Level Biomarker 1", "Age"), n_tot = c(0L, 0L), n_tot_events = c(0L, 0L), median = c(NA, NA), diff --git a/tests/testthat/test-prop_diff.R b/tests/testthat/test-prop_diff.R index fa3cd75ea5..fe6bf150f1 100644 --- a/tests/testthat/test-prop_diff.R +++ b/tests/testthat/test-prop_diff.R @@ -8,7 +8,7 @@ testthat::test_that("`prop_diff_ha` (proportion difference by Anderson-Hauck)", # according to SAS. expected <- list( diff = 0.25, - diff_ci = c(-0.9195, 1) + diff_ci = c(-0.9195, 1.0000) ) testthat::expect_equal(result, expected, tol = 0.0001) @@ -45,9 +45,8 @@ testthat::test_that("`prop_diff_nc` (proportion difference by Newcombe)", { # Edge case: Same proportion of response in A and B. rsp <- c(TRUE, FALSE, TRUE, FALSE) grp <- factor(c("A", "A", "B", "B"), levels = c("A", "B")) - result <- testthat::expect_warning( - prop_diff_nc(rsp = rsp, grp = grp, conf_level = 0.6) - ) + result <- prop_diff_nc(rsp = rsp, grp = grp, conf_level = 0.6) + # according to SAS. expected <- list( diff = 0, @@ -56,6 +55,59 @@ testthat::test_that("`prop_diff_nc` (proportion difference by Newcombe)", { testthat::expect_equal(result, expected, tol = 0.0001) }) +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")) + + result <- prop_diff_wald(rsp = rsp, grp = grp, conf_level = 0.9, correct = TRUE) + + expected <- list( + diff = 0.25, + diff_ci = c(-0.8069, 1.0000) + ) + testthat::expect_equal(result, expected, tol = 0.0001) + + # Edge case: Same proportion of response in A and B. + rsp <- c(TRUE, FALSE, TRUE, FALSE) + grp <- factor(c("A", "A", "B", "B"), levels = c("A", "B")) + result <- prop_diff_wald(rsp = rsp, grp = grp, conf_level = 0.6, correct = TRUE) + + expected <- list( + diff = 0, + diff_ci = c(-0.9208, 0.9208) + ) + testthat::expect_equal(result, expected, tol = 0.0001) +}) + +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")) + + result <- suppressWarnings( + prop_diff_wald(rsp = rsp, grp = grp, conf_level = 0.9, correct = FALSE) + ) + expected <- list( + diff = 0.25, + diff_ci = c(-0.4319, 0.9319) + ) + testthat::expect_equal(result, expected, tol = 0.0001) + + # Edge case: Same proportion of response in A and B. + rsp <- c(TRUE, FALSE, TRUE, FALSE) + grp <- factor(c("A", "A", "B", "B"), levels = c("A", "B")) + result <- prop_diff_wald(rsp = rsp, grp = grp, conf_level = 0.6, correct = FALSE) + + expected <- list( + diff = 0, + diff_ci = c(-0.4208, 0.4208) + ) + testthat::expect_equal(result, expected, tol = 0.0001) +}) + testthat::test_that("`prop_diff_cmh` (proportion difference by CMH)", { set.seed(2, kind = "Mersenne-Twister") diff --git a/tests/testthat/test-response_biomarkers_subgroups.R b/tests/testthat/test-response_biomarkers_subgroups.R index ebc89b76cc..fe8b58498b 100644 --- a/tests/testthat/test-response_biomarkers_subgroups.R +++ b/tests/testthat/test-response_biomarkers_subgroups.R @@ -1,5 +1,6 @@ library(scda) library(dplyr) +library(survival) preprocess_adrs <- function(adrs) { @@ -17,7 +18,7 @@ preprocess_adrs <- function(adrs) { ) } -adrs <- synthetic_cdisc_data("rcd_2021_05_05")$adrs +adrs <- synthetic_cdisc_data("rcd_2022_02_28")$adrs # extract_rsp_biomarkers ---- @@ -39,10 +40,10 @@ testthat::test_that("extract_rsp_biomarkers functions as expected with valid inp "AGE", "BMRKR1", "AGE", "BMRKR1", "AGE", "BMRKR1", "AGE", "BMRKR1" ), biomarker_label = c( - "Age", "Continous Level Biomarker 1", - "Age", "Continous Level Biomarker 1", "Age", "Continous Level Biomarker 1", - "Age", "Continous Level Biomarker 1", "Age", "Continous Level Biomarker 1", - "Age", "Continous Level Biomarker 1" + "Age", "Continuous Level Biomarker 1", + "Age", "Continuous Level Biomarker 1", "Age", "Continuous Level Biomarker 1", + "Age", "Continuous Level Biomarker 1", "Age", "Continuous Level Biomarker 1", + "Age", "Continuous Level Biomarker 1" ), n_tot = c( 400L, 400L, 231L, 231L, 169L, 169L, 135L, 135L, 135L, 135L, @@ -185,7 +186,7 @@ testthat::test_that("tabulate_rsp_biomarkers works as expected with valid input" ncol = 7, data = c( "", "Age", "All Patients", "Sex", "F", "M", "Categorical Level Biomarker 2", - "LOW", "MEDIUM", "HIGH", "Continous Level Biomarker 1", "All Patients", + "LOW", "MEDIUM", "HIGH", "Continuous Level Biomarker 1", "All Patients", "Sex", "F", "M", "Categorical Level Biomarker 2", "LOW", "MEDIUM", "HIGH", "Total n", "", "400", "", "231", "169", "", "135", "135", "130", "", "400", "", "231", "169", "", "135", "135", "130", @@ -229,14 +230,14 @@ testthat::test_that("tabulate_rsp_biomarkers functions as expected with NULL sub result <- tabulate_rsp_biomarkers(df) result_matrix <- to_string_matrix(result) - expected_first_col <- c("", "Age", "All Patients", "Continous Level Biomarker 1", "All Patients") + expected_first_col <- c("", "Age", "All Patients", "Continuous Level Biomarker 1", "All Patients") testthat::expect_identical(result_matrix[, 1L], expected_first_col) }) testthat::test_that("tabulate_rsp_biomarkers works with only a single biomarker in the data frame", { df1 <- data.frame( biomarker = "BMRKR1", - biomarker_label = "Continous Level Biomarker 1", + biomarker_label = "Continuous Level Biomarker 1", n_tot = 400L, n_rsp = 282L, prop = 0.705, @@ -255,7 +256,7 @@ testthat::test_that("tabulate_rsp_biomarkers works with only a single biomarker result_matrix <- to_string_matrix(result) expected_matrix <- matrix( data = c( - "", "Continous Level Biomarker 1", "All Patients", + "", "Continuous Level Biomarker 1", "All Patients", "Total n", "", "400", "Responders", "", "282", "Response (%)", "", "70.5%", "Odds Ratio", "", "0.98", "95% CI", "", "(0.95, 1.01)", "p-value (Wald)", "", "0.3000" diff --git a/tests/testthat/test-stat.R b/tests/testthat/test-stat.R index 5e0c345f2e..fa82857c7a 100644 --- a/tests/testthat/test-stat.R +++ b/tests/testthat/test-stat.R @@ -20,7 +20,8 @@ testthat::test_that("stat_mean_ci works for series without NAs testthat::expect_identical(result, expected) }) -testthat::test_that("stat_mean_ci works for series with NAs (including extreme case n = 1 and various n_min values)", { +testthat::test_that("stat_mean_ci works for series with NAs + (including extreme case n = 1 and various n_min values)", { # n = 0, na.rm = TRUE, n_min = 2 result <- stat_mean_ci(x = rep(NA, 10), gg_helper = FALSE) @@ -77,6 +78,89 @@ testthat::test_that("stat_mean_ci works for series with NAs (including extreme c testthat::expect_identical(result, expected) }) +testthat::test_that("stat_mean_pval works for series without NAs + (including extreme case n = 1 and various n_min values)", { + + # n = 1, na.rm = TRUE, n_min = 2 + result <- stat_mean_pval(x = 1) + expected <- c(p_value = NA_real_) + testthat::expect_identical(result, expected) + + # n = 2, na.rm = TRUE, n_min = 2 + result <- round(stat_mean_pval(x = 1:2), digits = 2) + expected <- c(p_value = 0.2) + testthat::expect_identical(result, expected) + + # n = 2, na.rm = TRUE, n_min = 3 + result <- stat_mean_pval(x = 1:2, n_min = 3) + expected <- c(p_value = NA_real_) + testthat::expect_identical(result, expected) +}) + +testthat::test_that("stat_mean_pval works for series with NAs + (including extreme case n = 1 and various n_min values)", { + + # n = 0, na.rm = TRUE, n_min = 2 + result <- stat_mean_pval(x = rep(NA, 10)) + expected <- c(p_value = NA_real_) + testthat::expect_identical(result, expected) + + # n = 0, na.rm = FALSE, n_min = 2 + result <- stat_mean_pval(x = rep(NA, 10), na.rm = FALSE) + expected <- c(p_value = NA_real_) + testthat::expect_identical(result, expected) + + # n = 1, na.rm = TRUE, n_min = 2 + result <- stat_mean_pval(x = c(1, NA, NA, NA)) + expected <- c(p_value = NA_real_) + testthat::expect_identical(result, expected) + + # n = 1, na.rm = FALSE, n_min = 2 + result <- stat_mean_pval(x = c(1, NA, NA, NA), na.rm = FALSE) + expected <- c(p_value = NA_real_) + testthat::expect_identical(result, expected) + + # n = 2, na.rm = TRUE, n_min = 2 + result <- round(stat_mean_pval(x = c(1:2, NA, NA, NA)), digits = 2) + expected <- c(p_value = 0.2) + testthat::expect_identical(result, expected) + + # n = 2, na.rm = TRUE, n_min = 3 + result <- stat_mean_pval(x = c(1:2, NA, NA, NA), n_min = 3) + expected <- c(p_value = NA_real_) + testthat::expect_identical(result, expected) + + # n = 2, na.rm = FALSE, n_min = 2 + result <- stat_mean_pval(x = c(1:2, NA, NA, NA), na.rm = FALSE) + expected <- c(p_value = NA_real_) + testthat::expect_identical(result, expected) + + # n = 2, na.rm = FALSE, n_min = 3 + result <- round( + stat_mean_pval( + x = c(1:2, NA, NA, NA), na.rm = FALSE, n_min = 3 + ), + digits = 2 + ) + expected <- c(p_value = NA_real_) + testthat::expect_identical(result, expected) +}) + +testthat::test_that("stat_mean_pval returns the correct p-value", { + + # n = 5, na.rm = TRUE, test_mean = 0 + x <- 1:5 + result <- stat_mean_pval(x) + expected <- c(p_value = t.test(x)$p.value) + testthat::expect_equal(result, expected) + + # n = 3, na.rm = TRUE, test_mean = 0.5 + x <- 1:3 + result <- stat_mean_pval(x, test_mean = 0.5) + expected <- c(p_value = t.test(x, mu = 0.5)$p.value) + testthat::expect_equal(result, expected) +}) + testthat::test_that("stat_median_ci works for series without NAs (including extreme case n = 1)", { # n = 1, na.rm = TRUE diff --git a/tests/testthat/test-summarize_change.R b/tests/testthat/test-summarize_change.R index 628920b1b2..b37ecc45fc 100644 --- a/tests/testthat/test-summarize_change.R +++ b/tests/testthat/test-summarize_change.R @@ -18,9 +18,11 @@ testthat::test_that("s_change_from_baseline handles empty data (complete missing sd = c(sd = NA_real_), se = c(se = NA_real_), mean_sd = c(mean = NA_real_, sd = NA_real_), + mean_se = c(mean = NA_real_, se = NA_real_), mean_ci = formatters::with_label(c(mean_ci_lwr = NA_real_, mean_ci_upr = NA_real_), "Mean 95% CI"), mean_sei = formatters::with_label(c(mean_sei_lwr = NA_real_, mean_sei_upr = NA_real_), "Mean -/+ 1xSE"), mean_sdi = formatters::with_label(c(mean_sdi_lwr = NA_real_, mean_sdi_upr = NA_real_), "Mean -/+ 1xSD"), + mean_pval = formatters::with_label(c(p_value = NA_real_), "Mean p-value (H0: mean = 0)"), median = c(median = NA_real_), mad = c(mad = NA_real_), median_ci = formatters::with_label( @@ -61,9 +63,11 @@ testthat::test_that("s_change_from_baseline handles NA in baseline values", { sd = c(sd = 3), se = c(se = 1.732051), mean_sd = c(mean = 3, sd = 3), + mean_se = c(mean = 3, se = 1.732051), mean_ci = formatters::with_label(c(mean_ci_lwr = -4.452413, mean_ci_upr = 10.452413), "Mean 95% CI"), mean_sei = formatters::with_label(c(mean_sei_lwr = 1.267949, mean_sei_upr = 4.732051), "Mean -/+ 1xSE"), mean_sdi = formatters::with_label(c(mean_sdi_lwr = 0, mean_sdi_upr = 6), "Mean -/+ 1xSD"), + mean_pval = formatters::with_label(c(p_value = 0.2254033), "Mean p-value (H0: mean = 0)"), median = c(median = 3), mad = c(mad = 0), median_ci = formatters::with_label( @@ -107,9 +111,11 @@ testthat::test_that("s_change_from_baseline handles baseline substitution", { sd = c(sd = 0.7071068), se = c(se = 0.5), mean_sd = c(mean = 1.5, sd = 0.7071068), + mean_se = c(mean = 1.5, se = 0.5), mean_ci = formatters::with_label(c(mean_ci_lwr = -4.853102, mean_ci_upr = 7.853102), "Mean 95% CI"), mean_sei = formatters::with_label(c(mean_sei_lwr = 1, mean_sei_upr = 2), "Mean -/+ 1xSE"), mean_sdi = formatters::with_label(c(mean_sdi_lwr = 0.7928932, mean_sdi_upr = 2.2071068), "Mean -/+ 1xSD"), + mean_pval = formatters::with_label(c(p_value = 0.2048328), "Mean p-value (H0: mean = 0)"), median = c(median = 1.5), mad = c(mad = 0), median_ci = formatters::with_label( @@ -137,9 +143,11 @@ testthat::test_that("s_change_from_baseline handles baseline substitution", { sd = c(sd = 2.12132), se = c(se = 1.5), mean_sd = c(mean = 2.5, sd = 2.12132), + mean_se = c(mean = 2.5, se = 1.5), mean_ci = formatters::with_label(c(mean_ci_lwr = -16.55931, mean_ci_upr = 21.55931), "Mean 95% CI"), mean_sei = formatters::with_label(c(mean_sei_lwr = 1, mean_sei_upr = 4), "Mean -/+ 1xSE"), mean_sdi = formatters::with_label(c(mean_sdi_lwr = 0.3786797, mean_sdi_upr = 4.6213203), "Mean -/+ 1xSD"), + mean_pval = formatters::with_label(c(p_value = 0.3440417), "Mean p-value (H0: mean = 0)"), median = c(median = 2.5), mad = c(mad = 0), median_ci = formatters::with_label( diff --git a/tests/testthat/test-summarize_variables.R b/tests/testthat/test-summarize_variables.R index 1c4aee9e44..ae046e2470 100644 --- a/tests/testthat/test-summarize_variables.R +++ b/tests/testthat/test-summarize_variables.R @@ -8,7 +8,8 @@ testthat::test_that("control_summarize_vars works with customized parameters", { expected <- list( conf_level = 0.9, quantiles = c(0.1, 0.9), - quantile_type = 2 + quantile_type = 2, + test_mean = 0 ) testthat::expect_identical(result, expected) }) @@ -29,9 +30,11 @@ testthat::test_that("s_summary return NA for x length 0L", { sd = c(sd = NA_real_), se = c(se = NA_real_), mean_sd = c(mean = NA_real_, sd = NA_real_), + mean_se = c(mean = NA_real_, se = NA_real_), mean_ci = formatters::with_label(c(mean_ci_lwr = NA_real_, mean_ci_upr = NA_real_), "Mean 95% CI"), mean_sei = formatters::with_label(c(mean_sei_lwr = NA_real_, mean_sei_upr = NA_real_), "Mean -/+ 1xSE"), mean_sdi = formatters::with_label(c(mean_sdi_lwr = NA_real_, mean_sdi_upr = NA_real_), "Mean -/+ 1xSD"), + mean_pval = formatters::with_label(c(p_value = NA_real_), "Mean p-value (H0: mean = 0)"), median = c(median = NA_real_), mad = c(mad = NA_real_), median_ci = formatters::with_label(c(median_ci_lwr = NA_real_, median_ci_upr = NA_real_), "Median 95% CI"), @@ -60,9 +63,11 @@ testthat::test_that("s_summary handles NA", { sd = c(sd = NA_real_), se = c(se = NA_real_), mean_sd = c(mean = 1, sd = NA_real_), + mean_se = c(mean = 1, se = NA_real_), mean_ci = formatters::with_label(c(mean_ci_lwr = NA_real_, mean_ci_upr = NA_real_), "Mean 95% CI"), mean_sei = formatters::with_label(c(mean_sei_lwr = NA_real_, mean_sei_upr = NA_real_), "Mean -/+ 1xSE"), mean_sdi = formatters::with_label(c(mean_sdi_lwr = NA_real_, mean_sdi_upr = NA_real_), "Mean -/+ 1xSD"), + mean_pval = formatters::with_label(c(p_value = NA_real_), "Mean p-value (H0: mean = 0)"), median = c(median = 1), mad = c(mad = 0), median_ci = formatters::with_label(c(median_ci_lwr = NA_real_, median_ci_upr = NA_real_), "Median 95% CI"), @@ -87,9 +92,11 @@ testthat::test_that("s_summary handles NA", { sd = c(sd = NA_real_), se = c(se = NA_real_), mean_sd = c(mean = NA_real_, sd = NA_real_), + mean_se = c(mean = NA_real_, se = NA_real_), mean_ci = formatters::with_label(c(mean_ci_lwr = NA_real_, mean_ci_upr = NA_real_), "Mean 95% CI"), mean_sei = formatters::with_label(c(mean_sei_lwr = NA_real_, mean_sei_upr = NA_real_), "Mean -/+ 1xSE"), mean_sdi = formatters::with_label(c(mean_sdi_lwr = NA_real_, mean_sdi_upr = NA_real_), "Mean -/+ 1xSD"), + mean_pval = formatters::with_label(c(p_value = NA_real_), "Mean p-value (H0: mean = 0)"), median = c(median = NA_real_), mad = c(mad = NA_real_), median_ci = formatters::with_label(c(median_ci_lwr = NA_real_, median_ci_upr = NA_real_), "Median 95% CI"), @@ -116,9 +123,11 @@ testthat::test_that("s_summary returns right results for n = 2", { sd = c(sd = 0.7071068), se = c(se = 0.5), mean_sd = c(mean = 1.5, sd = 0.7071068), + mean_se = c(mean = 1.5, se = 0.5), mean_ci = formatters::with_label(c(mean_ci_lwr = -4.853102, mean_ci_upr = 7.853102), "Mean 95% CI"), mean_sei = formatters::with_label(c(mean_sei_lwr = 1, mean_sei_upr = 2), "Mean -/+ 1xSE"), mean_sdi = formatters::with_label(c(mean_sdi_lwr = 0.7928932, mean_sdi_upr = 2.2071068), "Mean -/+ 1xSD"), + mean_pval = formatters::with_label(c(p_value = 0.2048328), "Mean p-value (H0: mean = 0)"), median = c(median = 1.5), mad = c(mad = 0), median_ci = formatters::with_label(c(median_ci_lwr = NA_real_, median_ci_upr = NA_real_), "Median 95% CI"), @@ -148,9 +157,11 @@ testthat::test_that("s_summary returns right results for n = 8", { sd = c(sd = 3.207135), se = c(se = 1.133893), mean_sd = c(mean = 6, sd = 3.207135), + mean_se = c(mean = 6, se = 1.133893), mean_ci = formatters::with_label(c(mean_ci_lwr = 3.318768, mean_ci_upr = 8.681232), "Mean 95% CI"), mean_sei = formatters::with_label(c(mean_sei_lwr = 4.866107, mean_sei_upr = 7.133893), "Mean -/+ 1xSE"), mean_sdi = formatters::with_label(c(mean_sdi_lwr = 2.792865, mean_sdi_upr = 9.207135), "Mean -/+ 1xSD"), + mean_pval = formatters::with_label(c(p_value = 0.001133783), "Mean p-value (H0: mean = 0)"), median = c(median = 6.5), mad = c(mad = 0), median_ci = formatters::with_label(c(median_ci_lwr = 1, median_ci_upr = 10), "Median 95% CI"), @@ -434,16 +445,16 @@ testthat::test_that("`summarize_vars` works with healthy input, and control func summarize_vars( vars = "AVAL", control = control_summarize_vars(quantiles = c(0.1, 0.9), conf_level = 0.9), - .stats = c("n", "mean_sd", "mean_ci", "quantiles") + .stats = c("n", "mean_sd", "mean_se", "mean_ci", "quantiles") ) result <- build_table(l, df = dta_test) expected <- structure( c( - "", "n", "Mean (SD)", "Mean 90% CI", "10% and 90%-ile", "all obs", - "9", "5.0 (2.7)", "(3.30, 6.70)", "1.0 - 9.0" + "", "n", "Mean (SD)", "Mean (SE)", "Mean 90% CI", "10% and 90%-ile", "all obs", + "9", "5.0 (2.7)", "5.0 (0.9)", "(3.30, 6.70)", "1.0 - 9.0" ), - .Dim = c(5L, 2L) + .Dim = c(6L, 2L) ) testthat::expect_identical(to_string_matrix(result), expected) diff --git a/tests/testthat/test-survival_biomarkers_subgroups.R b/tests/testthat/test-survival_biomarkers_subgroups.R index 1184821266..9f4ad72b16 100644 --- a/tests/testthat/test-survival_biomarkers_subgroups.R +++ b/tests/testthat/test-survival_biomarkers_subgroups.R @@ -20,7 +20,7 @@ preprocess_adtte <- function(adtte) { ) } -adtte <- synthetic_cdisc_data("rcd_2021_05_05")$adtte +adtte <- synthetic_cdisc_data("rcd_2022_02_28")$adtte # extract_survival_biomarkers ---- @@ -42,10 +42,10 @@ testthat::test_that("extract_survival_biomarkers functions as expected with vali "AGE", "BMRKR1", "AGE", "BMRKR1", "AGE", "BMRKR1", "AGE", "BMRKR1", "AGE", "BMRKR1", "AGE", "BMRKR1" ), biomarker_label = c( - "Age", "Continous Level Biomarker 1", - "Age", "Continous Level Biomarker 1", "Age", "Continous Level Biomarker 1", - "Age", "Continous Level Biomarker 1", "Age", "Continous Level Biomarker 1", - "Age", "Continous Level Biomarker 1" + "Age", "Continuous Level Biomarker 1", + "Age", "Continuous Level Biomarker 1", "Age", "Continuous Level Biomarker 1", + "Age", "Continuous Level Biomarker 1", "Age", "Continuous Level Biomarker 1", + "Age", "Continuous Level Biomarker 1" ), n_tot = c( 400L, 400L, 231L, 231L, 169L, 169L, 135L, 135L, 135L, 135L, 130L, 130L @@ -176,7 +176,7 @@ testthat::test_that("tabulate_survival_biomarkers works as expected with valid i ncol = 7, data = c( "", "Age", "All Patients", "Sex", "F", "M", "Categorical Level Biomarker 2", - "LOW", "MEDIUM", "HIGH", "Continous Level Biomarker 1", "All Patients", + "LOW", "MEDIUM", "HIGH", "Continuous Level Biomarker 1", "All Patients", "Sex", "F", "M", "Categorical Level Biomarker 2", "LOW", "MEDIUM", "HIGH", "Total n", "", "400", "", "231", "169", "", "135", "135", "130", "", "400", "", "231", "169", "", "135", "135", "130", @@ -221,14 +221,14 @@ testthat::test_that("tabulate_survival_biomarkers functions as expected with NUL result <- tabulate_survival_biomarkers(df, time_unit = as.character(adtte_f$AVALU[1])) result_matrix <- to_string_matrix(result) - expected_first_col <- c("", "Age", "All Patients", "Continous Level Biomarker 1", "All Patients") + expected_first_col <- c("", "Age", "All Patients", "Continuous Level Biomarker 1", "All Patients") testthat::expect_identical(result_matrix[, 1L], expected_first_col) }) testthat::test_that("tabulate_survival_biomarkers works with only a single biomarker in the data frame", { df1 <- data.frame( biomarker = "BMRKR1", - biomarker_label = "Continous Level Biomarker 1", + biomarker_label = "Continuous Level Biomarker 1", n_tot = 400L, n_tot_events = 282, median = 680, @@ -247,7 +247,7 @@ testthat::test_that("tabulate_survival_biomarkers works with only a single bioma result_matrix <- to_string_matrix(result) expected_matrix <- matrix( data = c( - "", "Continous Level Biomarker 1", "All Patients", + "", "Continuous Level Biomarker 1", "All Patients", "Total n", "", "400", "Total Events", "", "282", "Median", "", "680.0", "Hazard Ratio", "", "0.98", "95% Wald CI", "", "(0.95, 1.01)", "p-value (Wald)", "", "0.3000" diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R index e788d35729..e18b36913b 100644 --- a/tests/testthat/test-utils.R +++ b/tests/testthat/test-utils.R @@ -3,11 +3,23 @@ testthat::test_that("f_conf_level works for proportion", { expected <- "95% CI" testthat::expect_identical(result, expected) }) + testthat::test_that("f_conf_level fails for non-proportion input", { testthat::expect_error(f_conf_level(1.1)) testthat::expect_error(f_conf_level(-1)) }) +testthat::test_that("f_pval works for numeric input", { + result <- f_pval(0.05) + expected <- "p-value (H0: mean = 0.05)" + testthat::expect_identical(result, expected) +}) + +testthat::test_that("f_pval fails for non-numeric input", { + testthat::expect_error(f_pval("a")) + testthat::expect_error(f_pval(NULL)) +}) + testthat::test_that("make_names works as expected", { nams <- c("Any Grade (%)", "Total AE numbers!", "No adverse events ...") result <- make_names(nams)