Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Refactor estimate_incidence_rate #1300

Merged
merged 16 commits into from
Sep 20, 2024
Merged
Show file tree
Hide file tree
Changes from 15 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -130,6 +130,7 @@ Collate:
'h_adsl_adlb_merge_using_worst_flag.R'
'h_biomarkers_subgroups.R'
'h_cox_regression.R'
'h_incidence_rate.R'
'h_km.R'
'h_logistic_regression.R'
'h_map_for_count_abnormal.R'
Expand Down
1 change: 1 addition & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@ export(a_count_patients_with_event)
export(a_count_patients_with_flags)
export(a_count_values)
export(a_coxreg)
export(a_incidence_rate)
export(a_length_proportion)
export(a_odds_ratio)
export(a_proportion)
Expand Down
3 changes: 3 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,8 @@
* Added the `riskdiff` argument to `tabulate_rsp_subgroups` and `tabulate_survival_subgroups` to allow users to add a risk difference table column, and function `control_riskdiff` to specify settings for the risk difference column.
* Added warning to `tabulate_rsp_subgroups` when `pval` statistic is selected but `df` has not been correctly generated to add p-values to the output table.
* Added `n_rate` statistic as a non-default option to `estimate_incidence_rate` which returns both number of events observed and estimated incidence rate.
* Added `n_unique` statistic as a non-default option to `estimate_incidence_rate` which returns total number of patients with at least one event observed.
* Refactored `estimate_incidence_rate` to work as both an analyze function and a summarize function, controlled by the added `summarize` parameter. When `summarize = TRUE`, labels can be fine-tuned via the new `label_fmt` argument to the same function.

### Bug Fixes
* Fixed a bug in `a_surv_time` that threw an error when split only has `"is_event"`.
Expand All @@ -22,6 +24,7 @@
* Began deprecation of the confusing functions `summary_formats` and `summary_labels`.
* Enhanced general descriptions of analyze and summarize functions throughout package documentation.
* Finalized deprecation of the `strata` and `cohort_id` arguments to `g_lineplot`.
* Moved incidence rate helper functions into a separate `h_incidence_rate.R` file.

# tern 0.9.5

Expand Down
135 changes: 135 additions & 0 deletions R/h_incidence_rate.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,135 @@
#' Helper functions for incidence rate
#'
#' @description `r lifecycle::badge("stable")`
#'
#' @param control (`list`)\cr parameters for estimation details, specified by using
#' the helper function [control_incidence_rate()]. Possible parameter options are:
#' * `conf_level`: (`proportion`)\cr confidence level for the estimated incidence rate.
#' * `conf_type`: (`string`)\cr `normal` (default), `normal_log`, `exact`, or `byar`
#' for confidence interval type.
#' * `input_time_unit`: (`string`)\cr `day`, `week`, `month`, or `year` (default)
#' indicating time unit for data input.
#' * `num_pt_year`: (`numeric`)\cr time unit for desired output (in person-years).
#' @param person_years (`numeric(1)`)\cr total person-years at risk.
#' @param alpha (`numeric(1)`)\cr two-sided alpha-level for confidence interval.
#' @param n_events (`integer(1)`)\cr number of events observed.
#'
#' @return Estimated incidence rate, `rate`, and associated confidence interval, `rate_ci`.
#'
#' @seealso [incidence_rate]
#'
#' @name h_incidence_rate
NULL

#' @describeIn h_incidence_rate Helper function to estimate the incidence rate and
#' associated confidence interval.
#'
#' @keywords internal
h_incidence_rate <- function(person_years,
n_events,
control = control_incidence_rate()) {
alpha <- 1 - control$conf_level
est <- switch(control$conf_type,
normal = h_incidence_rate_normal(person_years, n_events, alpha),
normal_log = h_incidence_rate_normal_log(person_years, n_events, alpha),
exact = h_incidence_rate_exact(person_years, n_events, alpha),
byar = h_incidence_rate_byar(person_years, n_events, alpha)
)

num_pt_year <- control$num_pt_year
list(
rate = est$rate * num_pt_year,
rate_ci = est$rate_ci * num_pt_year
)
}

#' @describeIn h_incidence_rate Helper function to estimate the incidence rate and
#' associated confidence interval based on the normal approximation for the
#' incidence rate. Unit is one person-year.
#'
#' @examples
#' h_incidence_rate_normal(200, 2)
#'
#' @export
h_incidence_rate_normal <- function(person_years,
n_events,
alpha = 0.05) {
checkmate::assert_number(person_years)
checkmate::assert_number(n_events)
assert_proportion_value(alpha)

est <- n_events / person_years
se <- sqrt(est / person_years)
ci <- est + c(-1, 1) * stats::qnorm(1 - alpha / 2) * se

list(rate = est, rate_ci = ci)
}

#' @describeIn h_incidence_rate Helper function to estimate the incidence rate and
#' associated confidence interval based on the normal approximation for the
#' logarithm of the incidence rate. Unit is one person-year.
#'
#' @examples
#' h_incidence_rate_normal_log(200, 2)
#'
#' @export
h_incidence_rate_normal_log <- function(person_years,
n_events,
alpha = 0.05) {
checkmate::assert_number(person_years)
checkmate::assert_number(n_events)
assert_proportion_value(alpha)

rate_est <- n_events / person_years
rate_se <- sqrt(rate_est / person_years)
lrate_est <- log(rate_est)
lrate_se <- rate_se / rate_est
ci <- exp(lrate_est + c(-1, 1) * stats::qnorm(1 - alpha / 2) * lrate_se)

list(rate = rate_est, rate_ci = ci)
}

#' @describeIn h_incidence_rate Helper function to estimate the incidence rate and
#' associated exact confidence interval. Unit is one person-year.
#'
#' @examples
#' h_incidence_rate_exact(200, 2)
#'
#' @export
h_incidence_rate_exact <- function(person_years,
n_events,
alpha = 0.05) {
checkmate::assert_number(person_years)
checkmate::assert_number(n_events)
assert_proportion_value(alpha)

est <- n_events / person_years
lcl <- stats::qchisq(p = (alpha) / 2, df = 2 * n_events) / (2 * person_years)
ucl <- stats::qchisq(p = 1 - (alpha) / 2, df = 2 * n_events + 2) / (2 * person_years)

list(rate = est, rate_ci = c(lcl, ucl))
}

#' @describeIn h_incidence_rate Helper function to estimate the incidence rate and
#' associated Byar's confidence interval. Unit is one person-year.
#'
#' @examples
#' h_incidence_rate_byar(200, 2)
#'
#' @export
h_incidence_rate_byar <- function(person_years,
Melkiades marked this conversation as resolved.
Show resolved Hide resolved
n_events,
alpha = 0.05) {
checkmate::assert_number(person_years)
checkmate::assert_number(n_events)
assert_proportion_value(alpha)

est <- n_events / person_years
seg_1 <- n_events + 0.5
seg_2 <- 1 - 1 / (9 * (n_events + 0.5))
seg_3 <- stats::qnorm(1 - alpha / 2) * sqrt(1 / (n_events + 0.5)) / 3
lcl <- seg_1 * ((seg_2 - seg_3)^3) / person_years
edelarua marked this conversation as resolved.
Show resolved Hide resolved
ucl <- seg_1 * ((seg_2 + seg_3)^3) / person_years

list(rate = est, rate_ci = c(lcl, ucl))
}
Loading
Loading