1 |
#' Encode Categorical Missing Values in a Data Frame |
|
2 |
#' |
|
3 |
#' @description `r lifecycle::badge("stable")` |
|
4 |
#' |
|
5 |
#' This is a helper function to encode missing entries across groups of categorical |
|
6 |
#' variables in a data frame. |
|
7 |
#' |
|
8 |
#' @details Missing entries are those with `NA` or empty strings and will |
|
9 |
#' be replaced with a specified value. If factor variables include missing |
|
10 |
#' values, the missing value will be inserted as the last level. |
|
11 |
#' Similarly, in case character or logical variables should be converted to factors |
|
12 |
#' with the `char_as_factor` or `logical_as_factor` options, the missing values will |
|
13 |
#' be set as the last level. |
|
14 |
#' |
|
15 |
#' @param data (`data.frame`)\cr data set. |
|
16 |
#' @param omit_columns (`character`)\cr names of variables from `data` that should |
|
17 |
#' not be modified by this function. |
|
18 |
#' @param char_as_factor (`flag`)\cr whether to convert character variables |
|
19 |
#' in `data` to factors. |
|
20 |
#' @param logical_as_factor (`flag`)\cr whether to convert logical variables |
|
21 |
#' in `data` to factors. |
|
22 |
#' @param na_level (`string`)\cr used to replace all `NA` or empty |
|
23 |
#' values inside non-`omit_columns` columns. |
|
24 |
#' |
|
25 |
#' @return A `data.frame` with the chosen modifications applied. |
|
26 |
#' |
|
27 |
#' @seealso [sas_na()] and [explicit_na()] for other missing data helper functions. |
|
28 |
#' |
|
29 |
#' @examples |
|
30 |
#' my_data <- data.frame( |
|
31 |
#' u = c(TRUE, FALSE, NA, TRUE), |
|
32 |
#' v = factor(c("A", NA, NA, NA), levels = c("Z", "A")), |
|
33 |
#' w = c("A", "B", NA, "C"), |
|
34 |
#' x = c("D", "E", "F", NA), |
|
35 |
#' y = c("G", "H", "I", ""), |
|
36 |
#' z = c(1, 2, 3, 4), |
|
37 |
#' stringsAsFactors = FALSE |
|
38 |
#' ) |
|
39 |
#' |
|
40 |
#' # Example 1 |
|
41 |
#' # Encode missing values in all character or factor columns. |
|
42 |
#' df_explicit_na(my_data) |
|
43 |
#' # Also convert logical columns to factor columns. |
|
44 |
#' df_explicit_na(my_data, logical_as_factor = TRUE) |
|
45 |
#' # Encode missing values in a subset of columns. |
|
46 |
#' df_explicit_na(my_data, omit_columns = c("x", "y")) |
|
47 |
#' |
|
48 |
#' # Example 2 |
|
49 |
#' # Here we purposefully convert all `M` values to `NA` in the `SEX` variable. |
|
50 |
#' # After running `df_explicit_na` the `NA` values are encoded as `<Missing>` but they are not |
|
51 |
#' # included when generating `rtables`. |
|
52 |
#' adsl <- tern_ex_adsl |
|
53 |
#' adsl$SEX[adsl$SEX == "M"] <- NA |
|
54 |
#' adsl <- df_explicit_na(adsl) |
|
55 |
#' |
|
56 |
#' # If you want the `Na` values to be displayed in the table use the `na_level` argument. |
|
57 |
#' adsl <- tern_ex_adsl |
|
58 |
#' adsl$SEX[adsl$SEX == "M"] <- NA |
|
59 |
#' adsl <- df_explicit_na(adsl, na_level = "Missing Values") |
|
60 |
#' |
|
61 |
#' # Example 3 |
|
62 |
#' # Numeric variables that have missing values are not altered. This means that any `NA` value in |
|
63 |
#' # a numeric variable will not be included in the summary statistics, nor will they be included |
|
64 |
#' # in the denominator value for calculating the percent values. |
|
65 |
#' adsl <- tern_ex_adsl |
|
66 |
#' adsl$AGE[adsl$AGE < 30] <- NA |
|
67 |
#' adsl <- df_explicit_na(adsl) |
|
68 |
#' |
|
69 |
#' @export |
|
70 |
df_explicit_na <- function(data, |
|
71 |
omit_columns = NULL, |
|
72 |
char_as_factor = TRUE, |
|
73 |
logical_as_factor = FALSE, |
|
74 |
na_level = "<Missing>") { |
|
75 | 22x |
checkmate::assert_character(omit_columns, null.ok = TRUE, min.len = 1, any.missing = FALSE) |
76 | 21x |
checkmate::assert_data_frame(data) |
77 | 20x |
checkmate::assert_flag(char_as_factor) |
78 | 19x |
checkmate::assert_flag(logical_as_factor) |
79 | 19x |
checkmate::assert_string(na_level) |
80 | ||
81 | 17x |
target_vars <- if (is.null(omit_columns)) { |
82 | 15x |
names(data) |
83 |
} else { |
|
84 | 2x |
setdiff(names(data), omit_columns) # May have duplicates. |
85 |
} |
|
86 | 17x |
if (length(target_vars) == 0) { |
87 | 1x |
return(data) |
88 |
} |
|
89 | ||
90 | 16x |
l_target_vars <- split(target_vars, target_vars) |
91 | ||
92 |
# Makes sure target_vars exist in data and names are not duplicated. |
|
93 | 16x |
assert_df_with_variables(data, l_target_vars) |
94 | ||
95 | 16x |
for (x in target_vars) { |
96 | 304x |
xi <- data[[x]] |
97 | 304x |
xi_label <- obj_label(xi) |
98 | ||
99 |
# Determine whether to convert character or logical input. |
|
100 | 304x |
do_char_conversion <- is.character(xi) && char_as_factor |
101 | 304x |
do_logical_conversion <- is.logical(xi) && logical_as_factor |
102 | ||
103 |
# Pre-convert logical to character to deal correctly with replacing NA |
|
104 |
# values below. |
|
105 | 304x |
if (do_logical_conversion) { |
106 | 2x |
xi <- as.character(xi) |
107 |
} |
|
108 | ||
109 | 304x |
if (is.factor(xi) || is.character(xi)) { |
110 |
# Handle empty strings and NA values. |
|
111 | 217x |
xi <- explicit_na(sas_na(xi), label = na_level) |
112 | ||
113 |
# Convert to factors if requested for the original type, |
|
114 |
# set na_level as the last value. |
|
115 | 217x |
if (do_char_conversion || do_logical_conversion) { |
116 | 78x |
levels_xi <- setdiff(sort(unique(xi)), na_level) |
117 | 78x |
if (na_level %in% unique(xi)) { |
118 | 18x |
levels_xi <- c(levels_xi, na_level) |
119 |
} |
|
120 | ||
121 | 78x |
xi <- factor(xi, levels = levels_xi) |
122 |
} |
|
123 | ||
124 | 217x |
data[, x] <- formatters::with_label(xi, label = xi_label) |
125 |
} |
|
126 |
} |
|
127 | 16x |
return(data) |
128 |
} |
1 |
#' Odds Ratio Estimation |
|
2 |
#' |
|
3 |
#' @description `r lifecycle::badge("stable")` |
|
4 |
#' |
|
5 |
#' Compares bivariate responses between two groups in terms of odds ratios |
|
6 |
#' along with a confidence interval. |
|
7 |
#' |
|
8 |
#' @inheritParams argument_convention |
|
9 |
#' |
|
10 |
#' @details This function uses either logistic regression for unstratified |
|
11 |
#' analyses, or conditional logistic regression for stratified analyses. |
|
12 |
#' The Wald confidence interval with the specified confidence level is |
|
13 |
#' calculated. |
|
14 |
#' |
|
15 |
#' @note For stratified analyses, there is currently no implementation for conditional |
|
16 |
#' likelihood confidence intervals, therefore the likelihood confidence interval is not |
|
17 |
#' yet available as an option. Besides, when `rsp` contains only responders or non-responders, |
|
18 |
#' then the result values will be `NA`, because no odds ratio estimation is possible. |
|
19 |
#' |
|
20 |
#' @seealso Relevant helper function [h_odds_ratio()]. |
|
21 |
#' |
|
22 |
#' @name odds_ratio |
|
23 |
NULL |
|
24 | ||
25 |
#' @describeIn odds_ratio Statistics function which estimates the odds ratio |
|
26 |
#' between a treatment and a control. A `variables` list with `arm` and `strata` |
|
27 |
#' variable names must be passed if a stratified analysis is required. |
|
28 |
#' |
|
29 |
#' @inheritParams split_cols_by_groups |
|
30 |
#' |
|
31 |
#' @return |
|
32 |
#' * `s_odds_ratio()` returns a named list with the statistics `or_ci` |
|
33 |
#' (containing `est`, `lcl`, and `ucl`) and `n_tot`. |
|
34 |
#' |
|
35 |
#' @examples |
|
36 |
#' set.seed(12) |
|
37 |
#' dta <- data.frame( |
|
38 |
#' rsp = sample(c(TRUE, FALSE), 100, TRUE), |
|
39 |
#' grp = factor(rep(c("A", "B"), each = 50), levels = c("B", "A")), |
|
40 |
#' strata = factor(sample(c("C", "D"), 100, TRUE)) |
|
41 |
#' ) |
|
42 |
#' |
|
43 |
#' # Unstratified analysis. |
|
44 |
#' s_odds_ratio( |
|
45 |
#' df = subset(dta, grp == "A"), |
|
46 |
#' .var = "rsp", |
|
47 |
#' .ref_group = subset(dta, grp == "B"), |
|
48 |
#' .in_ref_col = FALSE, |
|
49 |
#' .df_row = dta |
|
50 |
#' ) |
|
51 |
#' |
|
52 |
#' # Stratified analysis. |
|
53 |
#' s_odds_ratio( |
|
54 |
#' df = subset(dta, grp == "A"), |
|
55 |
#' .var = "rsp", |
|
56 |
#' .ref_group = subset(dta, grp == "B"), |
|
57 |
#' .in_ref_col = FALSE, |
|
58 |
#' .df_row = dta, |
|
59 |
#' variables = list(arm = "grp", strata = "strata") |
|
60 |
#' ) |
|
61 |
#' |
|
62 |
#' @export |
|
63 |
s_odds_ratio <- function(df, |
|
64 |
.var, |
|
65 |
.ref_group, |
|
66 |
.in_ref_col, |
|
67 |
.df_row, |
|
68 |
variables = list(arm = NULL, strata = NULL), |
|
69 |
conf_level = 0.95, |
|
70 |
groups_list = NULL) { |
|
71 | 65x |
y <- list(or_ci = "", n_tot = "") |
72 | ||
73 | 65x |
if (!.in_ref_col) { |
74 | 65x |
assert_proportion_value(conf_level) |
75 | 65x |
assert_df_with_variables(df, list(rsp = .var)) |
76 | 65x |
assert_df_with_variables(.ref_group, list(rsp = .var)) |
77 | ||
78 | 65x |
if (is.null(variables$strata)) { |
79 | 52x |
data <- data.frame( |
80 | 52x |
rsp = c(.ref_group[[.var]], df[[.var]]), |
81 | 52x |
grp = factor( |
82 | 52x |
rep(c("ref", "Not-ref"), c(nrow(.ref_group), nrow(df))), |
83 | 52x |
levels = c("ref", "Not-ref") |
84 |
) |
|
85 |
) |
|
86 | 52x |
y <- or_glm(data, conf_level = conf_level) |
87 |
} else { |
|
88 | 13x |
assert_df_with_variables(.df_row, c(list(rsp = .var), variables)) |
89 | ||
90 |
# The group variable prepared for clogit must be synchronised with combination groups definition. |
|
91 | 13x |
if (is.null(groups_list)) { |
92 | 12x |
ref_grp <- as.character(unique(.ref_group[[variables$arm]])) |
93 | 12x |
trt_grp <- as.character(unique(df[[variables$arm]])) |
94 | 12x |
grp <- stats::relevel(factor(.df_row[[variables$arm]]), ref = ref_grp) |
95 |
} else { |
|
96 |
# If more than one level in reference col. |
|
97 | 1x |
reference <- as.character(unique(.ref_group[[variables$arm]])) |
98 | 1x |
grp_ref_flag <- vapply( |
99 | 1x |
X = groups_list, |
100 | 1x |
FUN.VALUE = TRUE, |
101 | 1x |
FUN = function(x) all(reference %in% x) |
102 |
) |
|
103 | 1x |
ref_grp <- names(groups_list)[grp_ref_flag] |
104 | ||
105 |
# If more than one level in treatment col. |
|
106 | 1x |
treatment <- as.character(unique(df[[variables$arm]])) |
107 | 1x |
grp_trt_flag <- vapply( |
108 | 1x |
X = groups_list, |
109 | 1x |
FUN.VALUE = TRUE, |
110 | 1x |
FUN = function(x) all(treatment %in% x) |
111 |
) |
|
112 | 1x |
trt_grp <- names(groups_list)[grp_trt_flag] |
113 | ||
114 | 1x |
grp <- combine_levels(.df_row[[variables$arm]], levels = reference, new_level = ref_grp) |
115 | 1x |
grp <- combine_levels(grp, levels = treatment, new_level = trt_grp) |
116 |
} |
|
117 | ||
118 |
# The reference level in `grp` must be the same as in the `rtables` column split. |
|
119 | 13x |
data <- data.frame( |
120 | 13x |
rsp = .df_row[[.var]], |
121 | 13x |
grp = grp, |
122 | 13x |
strata = interaction(.df_row[variables$strata]) |
123 |
) |
|
124 | 13x |
y_all <- or_clogit(data, conf_level = conf_level) |
125 | 13x |
checkmate::assert_string(trt_grp) |
126 | 13x |
checkmate::assert_subset(trt_grp, names(y_all$or_ci)) |
127 | 12x |
y$or_ci <- y_all$or_ci[[trt_grp]] |
128 | 12x |
y$n_tot <- y_all$n_tot |
129 |
} |
|
130 |
} |
|
131 | ||
132 | 64x |
y$or_ci <- formatters::with_label( |
133 | 64x |
x = y$or_ci, |
134 | 64x |
label = paste0("Odds Ratio (", 100 * conf_level, "% CI)") |
135 |
) |
|
136 | ||
137 | 64x |
y$n_tot <- formatters::with_label( |
138 | 64x |
x = y$n_tot, |
139 | 64x |
label = "Total n" |
140 |
) |
|
141 | ||
142 | 64x |
y |
143 |
} |
|
144 | ||
145 |
#' @describeIn odds_ratio Formatted analysis function which is used as `afun` in `estimate_odds_ratio()`. |
|
146 |
#' |
|
147 |
#' @return |
|
148 |
#' * `a_odds_ratio()` returns the corresponding list with formatted [rtables::CellValue()]. |
|
149 |
#' |
|
150 |
#' @examples |
|
151 |
#' a_odds_ratio( |
|
152 |
#' df = subset(dta, grp == "A"), |
|
153 |
#' .var = "rsp", |
|
154 |
#' .ref_group = subset(dta, grp == "B"), |
|
155 |
#' .in_ref_col = FALSE, |
|
156 |
#' .df_row = dta |
|
157 |
#' ) |
|
158 |
#' |
|
159 |
#' @export |
|
160 |
a_odds_ratio <- make_afun( |
|
161 |
s_odds_ratio, |
|
162 |
.formats = c(or_ci = "xx.xx (xx.xx - xx.xx)"), |
|
163 |
.indent_mods = c(or_ci = 1L) |
|
164 |
) |
|
165 | ||
166 |
#' @describeIn odds_ratio Layout-creating function which can take statistics function arguments |
|
167 |
#' and additional format arguments. This function is a wrapper for [rtables::analyze()]. |
|
168 |
#' |
|
169 |
#' @param ... arguments passed to `s_odds_ratio()`. |
|
170 |
#' |
|
171 |
#' @return |
|
172 |
#' * `estimate_odds_ratio()` returns a layout object suitable for passing to further layouting functions, |
|
173 |
#' or to [rtables::build_table()]. Adding this function to an `rtable` layout will add formatted rows containing |
|
174 |
#' the statistics from `s_odds_ratio()` to the table layout. |
|
175 |
#' |
|
176 |
#' @examples |
|
177 |
#' dta <- data.frame( |
|
178 |
#' rsp = sample(c(TRUE, FALSE), 100, TRUE), |
|
179 |
#' grp = factor(rep(c("A", "B"), each = 50)) |
|
180 |
#' ) |
|
181 |
#' |
|
182 |
#' l <- basic_table() %>% |
|
183 |
#' split_cols_by(var = "grp", ref_group = "B") %>% |
|
184 |
#' estimate_odds_ratio(vars = "rsp") |
|
185 |
#' |
|
186 |
#' build_table(l, df = dta) |
|
187 |
#' |
|
188 |
#' @export |
|
189 |
estimate_odds_ratio <- function(lyt, |
|
190 |
vars, |
|
191 |
nested = TRUE, |
|
192 |
..., |
|
193 |
show_labels = "hidden", |
|
194 |
table_names = vars, |
|
195 |
.stats = "or_ci", |
|
196 |
.formats = NULL, |
|
197 |
.labels = NULL, |
|
198 |
.indent_mods = NULL) { |
|
199 | 3x |
afun <- make_afun( |
200 | 3x |
a_odds_ratio, |
201 | 3x |
.stats = .stats, |
202 | 3x |
.formats = .formats, |
203 | 3x |
.labels = .labels, |
204 | 3x |
.indent_mods = .indent_mods |
205 |
) |
|
206 | ||
207 | 3x |
analyze( |
208 | 3x |
lyt, |
209 | 3x |
vars, |
210 | 3x |
afun = afun, |
211 | 3x |
nested = nested, |
212 | 3x |
extra_args = list(...), |
213 | 3x |
show_labels = show_labels, |
214 | 3x |
table_names = table_names |
215 |
) |
|
216 |
} |
|
217 | ||
218 |
#' Helper Functions for Odds Ratio Estimation |
|
219 |
#' |
|
220 |
#' @description `r lifecycle::badge("stable")` |
|
221 |
#' |
|
222 |
#' Functions to calculate odds ratios in [estimate_odds_ratio()]. |
|
223 |
#' |
|
224 |
#' @inheritParams argument_convention |
|
225 |
#' @param data (`data.frame`)\cr data frame containing at least the variables `rsp` and `grp`, and optionally |
|
226 |
#' `strata` for [or_clogit()]. |
|
227 |
#' |
|
228 |
#' @return A named `list` of elements `or_ci` and `n_tot`. |
|
229 |
#' |
|
230 |
#' @seealso [odds_ratio] |
|
231 |
#' |
|
232 |
#' @name h_odds_ratio |
|
233 |
NULL |
|
234 | ||
235 |
#' @describeIn h_odds_ratio Estimates the odds ratio based on [stats::glm()]. Note that there must be |
|
236 |
#' exactly 2 groups in `data` as specified by the `grp` variable. |
|
237 |
#' |
|
238 |
#' @examples |
|
239 |
#' # Data with 2 groups. |
|
240 |
#' data <- data.frame( |
|
241 |
#' rsp = as.logical(c(1, 1, 0, 1, 0, 0, 1, 1)), |
|
242 |
#' grp = letters[c(1, 1, 1, 2, 2, 2, 1, 2)], |
|
243 |
#' strata = letters[c(1, 2, 1, 2, 2, 2, 1, 2)], |
|
244 |
#' stringsAsFactors = TRUE |
|
245 |
#' ) |
|
246 |
#' |
|
247 |
#' # Odds ratio based on glm. |
|
248 |
#' or_glm(data, conf_level = 0.95) |
|
249 |
#' |
|
250 |
#' @export |
|
251 |
or_glm <- function(data, conf_level) { |
|
252 | 55x |
checkmate::assert_logical(data$rsp) |
253 | 55x |
assert_proportion_value(conf_level) |
254 | 55x |
assert_df_with_variables(data, list(rsp = "rsp", grp = "grp")) |
255 | 55x |
checkmate::assert_multi_class(data$grp, classes = c("factor", "character")) |
256 | ||
257 | 55x |
data$grp <- as_factor_keep_attributes(data$grp) |
258 | 55x |
assert_df_with_factors(data, list(val = "grp"), min.levels = 2, max.levels = 2) |
259 | 55x |
formula <- stats::as.formula("rsp ~ grp") |
260 | 55x |
model_fit <- stats::glm( |
261 | 55x |
formula = formula, data = data, |
262 | 55x |
family = stats::binomial(link = "logit") |
263 |
) |
|
264 | ||
265 |
# Note that here we need to discard the intercept. |
|
266 | 55x |
or <- exp(stats::coef(model_fit)[-1]) |
267 | 55x |
or_ci <- exp( |
268 | 55x |
stats::confint.default(model_fit, level = conf_level)[-1, , drop = FALSE] |
269 |
) |
|
270 | ||
271 | 55x |
values <- stats::setNames(c(or, or_ci), c("est", "lcl", "ucl")) |
272 | 55x |
n_tot <- stats::setNames(nrow(model_fit$model), "n_tot") |
273 | ||
274 | 55x |
list(or_ci = values, n_tot = n_tot) |
275 |
} |
|
276 | ||
277 |
#' @describeIn h_odds_ratio estimates the odds ratio based on [survival::clogit()]. This is done for |
|
278 |
#' the whole data set including all groups, since the results are not the same as when doing |
|
279 |
#' pairwise comparisons between the groups. |
|
280 |
#' |
|
281 |
#' @examples |
|
282 |
#' # Data with 3 groups. |
|
283 |
#' data <- data.frame( |
|
284 |
#' rsp = as.logical(c(1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0)), |
|
285 |
#' grp = letters[c(1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3)], |
|
286 |
#' strata = LETTERS[c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2)], |
|
287 |
#' stringsAsFactors = TRUE |
|
288 |
#' ) |
|
289 |
#' |
|
290 |
#' # Odds ratio based on stratified estimation by conditional logistic regression. |
|
291 |
#' or_clogit(data, conf_level = 0.95) |
|
292 |
#' |
|
293 |
#' @export |
|
294 |
or_clogit <- function(data, conf_level) { |
|
295 | 16x |
checkmate::assert_logical(data$rsp) |
296 | 16x |
assert_proportion_value(conf_level) |
297 | 16x |
assert_df_with_variables(data, list(rsp = "rsp", grp = "grp", strata = "strata")) |
298 | 16x |
checkmate::assert_multi_class(data$grp, classes = c("factor", "character")) |
299 | 16x |
checkmate::assert_multi_class(data$strata, classes = c("factor", "character")) |
300 | ||
301 | 16x |
data$grp <- as_factor_keep_attributes(data$grp) |
302 | 16x |
data$strata <- as_factor_keep_attributes(data$strata) |
303 | ||
304 |
# Deviation from convention: `survival::strata` must be simply `strata`. |
|
305 | 16x |
formula <- stats::as.formula("rsp ~ grp + strata(strata)") |
306 | 16x |
model_fit <- clogit_with_tryCatch(formula = formula, data = data) |
307 | ||
308 |
# Create a list with one set of OR estimates and CI per coefficient, i.e. |
|
309 |
# comparison of one group vs. the reference group. |
|
310 | 16x |
coef_est <- stats::coef(model_fit) |
311 | 16x |
ci_est <- stats::confint(model_fit, level = conf_level) |
312 | 16x |
or_ci <- list() |
313 | 16x |
for (coef_name in names(coef_est)) { |
314 | 18x |
grp_name <- gsub("^grp", "", x = coef_name) |
315 | 18x |
or_ci[[grp_name]] <- stats::setNames( |
316 | 18x |
object = exp(c(coef_est[coef_name], ci_est[coef_name, , drop = TRUE])), |
317 | 18x |
nm = c("est", "lcl", "ucl") |
318 |
) |
|
319 |
} |
|
320 | 16x |
list(or_ci = or_ci, n_tot = c(n_tot = model_fit$n)) |
321 |
} |
1 |
#' Incidence Rate |
|
2 |
#' |
|
3 |
#' @description `r lifecycle::badge("stable")` |
|
4 |
#' |
|
5 |
#' Estimate the event rate adjusted for person-years at risk, otherwise known |
|
6 |
#' as incidence rate. Primary analysis variable is the person-years at risk. |
|
7 |
#' |
|
8 |
#' @inheritParams argument_convention |
|
9 |
#' @param control (`list`)\cr parameters for estimation details, specified by using |
|
10 |
#' the helper function [control_incidence_rate()]. Possible parameter options are: |
|
11 |
#' * `conf_level` (`proportion`)\cr confidence level for the estimated incidence rate. |
|
12 |
#' * `conf_type` (`string`)\cr `normal` (default), `normal_log`, `exact`, or `byar` |
|
13 |
#' for confidence interval type. |
|
14 |
#' * `input_time_unit` (`string`)\cr `day`, `week`, `month`, or `year` (default) |
|
15 |
#' indicating time unit for data input. |
|
16 |
#' * `num_pt_year` (`numeric`)\cr time unit for desired output (in person-years). |
|
17 |
#' @param n_events (`integer`)\cr number of events observed. |
|
18 |
#' |
|
19 |
#' @seealso [control_incidence_rate()] and helper functions [h_incidence_rate]. |
|
20 |
#' |
|
21 |
#' @name incidence_rate |
|
22 |
NULL |
|
23 | ||
24 |
#' @describeIn incidence_rate Statistics function which estimates the incidence rate and the |
|
25 |
#' associated confidence interval. |
|
26 |
#' |
|
27 |
#' @return |
|
28 |
#' * `s_incidence_rate()` returns the following statistics: |
|
29 |
#' - `person_years`: Total person-years at risk. |
|
30 |
#' - `n_events`: Total number of events observed. |
|
31 |
#' - `rate`: Estimated incidence rate. |
|
32 |
#' - `rate_ci`: Confidence interval for the incidence rate. |
|
33 |
#' |
|
34 |
#' @examples |
|
35 |
#' library(dplyr) |
|
36 |
#' |
|
37 |
#' df <- data.frame( |
|
38 |
#' USUBJID = as.character(seq(6)), |
|
39 |
#' CNSR = c(0, 1, 1, 0, 0, 0), |
|
40 |
#' AVAL = c(10.1, 20.4, 15.3, 20.8, 18.7, 23.4), |
|
41 |
#' ARM = factor(c("A", "A", "A", "B", "B", "B")) |
|
42 |
#' ) %>% |
|
43 |
#' mutate(is_event = CNSR == 0) %>% |
|
44 |
#' mutate(n_events = as.integer(is_event)) |
|
45 |
#' |
|
46 |
#' @keywords internal |
|
47 |
s_incidence_rate <- function(df, |
|
48 |
.var, |
|
49 |
n_events, |
|
50 |
is_event, |
|
51 |
control = control_incidence_rate()) { |
|
52 | 1x |
if (!missing(is_event)) { |
53 | ! |
warning("argument is_event will be deprecated. Please use n_events.") |
54 | ||
55 | ! |
if (missing(n_events)) { |
56 | ! |
assert_df_with_variables(df, list(tte = .var, is_event = is_event)) |
57 | ! |
checkmate::assert_string(.var) |
58 | ! |
checkmate::assert_logical(df[[is_event]], any.missing = FALSE) |
59 | ! |
checkmate::assert_numeric(df[[.var]], any.missing = FALSE) |
60 | ! |
n_events <- is_event |
61 |
} |
|
62 |
} else { |
|
63 | 1x |
assert_df_with_variables(df, list(tte = .var, n_events = n_events)) |
64 | 1x |
checkmate::assert_string(.var) |
65 | 1x |
checkmate::assert_numeric(df[[.var]], any.missing = FALSE) |
66 | 1x |
checkmate::assert_integer(df[[n_events]], any.missing = FALSE) |
67 |
} |
|
68 | ||
69 | 1x |
input_time_unit <- control$input_time_unit |
70 | 1x |
num_pt_year <- control$num_pt_year |
71 | 1x |
conf_level <- control$conf_level |
72 | 1x |
person_years <- sum(df[[.var]], na.rm = TRUE) * ( |
73 | 1x |
1 * (input_time_unit == "year") + |
74 | 1x |
1 / 12 * (input_time_unit == "month") + |
75 | 1x |
1 / 52.14 * (input_time_unit == "week") + |
76 | 1x |
1 / 365.24 * (input_time_unit == "day") |
77 |
) |
|
78 | 1x |
n_events <- sum(df[[n_events]], na.rm = TRUE) |
79 | ||
80 | 1x |
result <- h_incidence_rate( |
81 | 1x |
person_years, |
82 | 1x |
n_events, |
83 | 1x |
control |
84 |
) |
|
85 | 1x |
list( |
86 | 1x |
person_years = formatters::with_label(person_years, "Total patient-years at risk"), |
87 | 1x |
n_events = formatters::with_label(n_events, "Number of adverse events observed"), |
88 | 1x |
rate = formatters::with_label(result$rate, paste("AE rate per", num_pt_year, "patient-years")), |
89 | 1x |
rate_ci = formatters::with_label(result$rate_ci, f_conf_level(conf_level)) |
90 |
) |
|
91 |
} |
|
92 | ||
93 |
#' @describeIn incidence_rate Formatted analysis function which is used as `afun` |
|
94 |
#' in `estimate_incidence_rate()`. |
|
95 |
#' |
|
96 |
#' @return |
|
97 |
#' * `a_incidence_rate()` returns the corresponding list with formatted [rtables::CellValue()]. |
|
98 |
#' |
|
99 |
#' |
|
100 |
#' @keywords internal |
|
101 |
a_incidence_rate <- make_afun( |
|
102 |
s_incidence_rate, |
|
103 |
.formats = c( |
|
104 |
"person_years" = "xx.x", |
|
105 |
"n_events" = "xx", |
|
106 |
"rate" = "xx.xx", |
|
107 |
"rate_ci" = "(xx.xx, xx.xx)" |
|
108 |
) |
|
109 |
) |
|
110 | ||
111 |
#' @describeIn incidence_rate Layout-creating function which can take statistics function arguments |
|
112 |
#' and additional format arguments. This function is a wrapper for [rtables::analyze()]. |
|
113 |
#' |
|
114 |
#' @return |
|
115 |
#' * `estimate_incidence_rate()` returns a layout object suitable for passing to further layouting functions, |
|
116 |
#' or to [rtables::build_table()]. Adding this function to an `rtable` layout will add formatted rows containing |
|
117 |
#' the statistics from `s_incidence_rate()` to the table layout. |
|
118 |
#' |
|
119 |
#' @examples |
|
120 |
#' basic_table() %>% |
|
121 |
#' split_cols_by("ARM") %>% |
|
122 |
#' add_colcounts() %>% |
|
123 |
#' estimate_incidence_rate( |
|
124 |
#' vars = "AVAL", |
|
125 |
#' n_events = "n_events", |
|
126 |
#' control = control_incidence_rate( |
|
127 |
#' input_time_unit = "month", |
|
128 |
#' num_pt_year = 100 |
|
129 |
#' ) |
|
130 |
#' ) %>% |
|
131 |
#' build_table(df) |
|
132 |
#' |
|
133 |
#' @export |
|
134 |
estimate_incidence_rate <- function(lyt, |
|
135 |
vars, |
|
136 |
nested = TRUE, |
|
137 |
..., |
|
138 |
show_labels = "hidden", |
|
139 |
table_names = vars, |
|
140 |
.stats = NULL, |
|
141 |
.formats = NULL, |
|
142 |
.labels = NULL, |
|
143 |
.indent_mods = NULL) { |
|
144 | 1x |
afun <- make_afun( |
145 | 1x |
a_incidence_rate, |
146 | 1x |
.stats = .stats, |
147 | 1x |
.formats = .formats, |
148 | 1x |
.labels = .labels, |
149 | 1x |
.indent_mods = .indent_mods |
150 |
) |
|
151 | ||
152 | 1x |
analyze( |
153 | 1x |
lyt, |
154 | 1x |
vars, |
155 | 1x |
show_labels = show_labels, |
156 | 1x |
table_names = table_names, |
157 | 1x |
afun = afun, |
158 | 1x |
nested = nested, |
159 | 1x |
extra_args = list(...) |
160 |
) |
|
161 |
} |
|
162 | ||
163 |
#' Helper Functions for Incidence Rate |
|
164 |
#' |
|
165 |
#' @description `r lifecycle::badge("stable")` |
|
166 |
#' |
|
167 |
#' @param control (`list`)\cr parameters for estimation details, specified by using |
|
168 |
#' the helper function [control_incidence_rate()]. Possible parameter options are: |
|
169 |
#' * `conf_level`: (`proportion`)\cr confidence level for the estimated incidence rate. |
|
170 |
#' * `conf_type`: (`string`)\cr `normal` (default), `normal_log`, `exact`, or `byar` |
|
171 |
#' for confidence interval type. |
|
172 |
#' * `input_time_unit`: (`string`)\cr `day`, `week`, `month`, or `year` (default) |
|
173 |
#' indicating time unit for data input. |
|
174 |
#' * `num_pt_year`: (`numeric`)\cr time unit for desired output (in person-years). |
|
175 |
#' @param person_years (`numeric`)\cr total person-years at risk. |
|
176 |
#' @param alpha (`numeric`)\cr two-sided alpha-level for confidence interval. |
|
177 |
#' @param n_events (`integer`)\cr number of events observed. |
|
178 |
#' |
|
179 |
#' @return Estimated incidence rate `rate` and associated confidence interval `rate_ci`. |
|
180 |
#' |
|
181 |
#' @seealso [incidence_rate] |
|
182 |
#' |
|
183 |
#' @name h_incidence_rate |
|
184 |
NULL |
|
185 | ||
186 |
#' @describeIn h_incidence_rate Helper function to estimate the incidence rate and |
|
187 |
#' associated confidence interval based on the normal approximation for the |
|
188 |
#' incidence rate. Unit is one person-year. |
|
189 |
#' |
|
190 |
#' @examples |
|
191 |
#' h_incidence_rate_normal(200, 2) |
|
192 |
#' |
|
193 |
#' @export |
|
194 |
h_incidence_rate_normal <- function(person_years, |
|
195 |
n_events, |
|
196 |
alpha = 0.05) { |
|
197 | 1x |
checkmate::assert_number(person_years) |
198 | 1x |
checkmate::assert_number(n_events) |
199 | 1x |
assert_proportion_value(alpha) |
200 | ||
201 | 1x |
est <- n_events / person_years |
202 | 1x |
se <- sqrt(est / person_years) |
203 | 1x |
ci <- est + c(-1, 1) * stats::qnorm(1 - alpha / 2) * se |
204 | ||
205 | 1x |
list(rate = est, rate_ci = ci) |
206 |
} |
|
207 | ||
208 |
#' @describeIn h_incidence_rate Helper function to estimate the incidence rate and |
|
209 |
#' associated confidence interval based on the normal approximation for the |
|
210 |
#' logarithm of the incidence rate. Unit is one person-year. |
|
211 |
#' |
|
212 |
#' @examples |
|
213 |
#' h_incidence_rate_normal_log(200, 2) |
|
214 |
#' |
|
215 |
#' @export |
|
216 |
h_incidence_rate_normal_log <- function(person_years, |
|
217 |
n_events, |
|
218 |
alpha = 0.05) { |
|
219 | 5x |
checkmate::assert_number(person_years) |
220 | 5x |
checkmate::assert_number(n_events) |
221 | 5x |
assert_proportion_value(alpha) |
222 | ||
223 | 5x |
rate_est <- n_events / person_years |
224 | 5x |
rate_se <- sqrt(rate_est / person_years) |
225 | 5x |
lrate_est <- log(rate_est) |
226 | 5x |
lrate_se <- rate_se / rate_est |
227 | 5x |
ci <- exp(lrate_est + c(-1, 1) * stats::qnorm(1 - alpha / 2) * lrate_se) |
228 | ||
229 | 5x |
list(rate = rate_est, rate_ci = ci) |
230 |
} |
|
231 | ||
232 |
#' @describeIn h_incidence_rate Helper function to estimate the incidence rate and |
|
233 |
#' associated exact confidence interval. Unit is one person-year. |
|
234 |
#' |
|
235 |
#' @examples |
|
236 |
#' h_incidence_rate_exact(200, 2) |
|
237 |
#' |
|
238 |
#' @export |
|
239 |
h_incidence_rate_exact <- function(person_years, |
|
240 |
n_events, |
|
241 |
alpha = 0.05) { |
|
242 | 1x |
checkmate::assert_number(person_years) |
243 | 1x |
checkmate::assert_number(n_events) |
244 | 1x |
assert_proportion_value(alpha) |
245 | ||
246 | 1x |
est <- n_events / person_years |
247 | 1x |
lcl <- stats::qchisq(p = (alpha) / 2, df = 2 * n_events) / (2 * person_years) |
248 | 1x |
ucl <- stats::qchisq(p = 1 - (alpha) / 2, df = 2 * n_events + 2) / (2 * person_years) |
249 | ||
250 | 1x |
list(rate = est, rate_ci = c(lcl, ucl)) |
251 |
} |
|
252 | ||
253 |
#' @describeIn h_incidence_rate Helper function to estimate the incidence rate and |
|
254 |
#' associated `Byar`'s confidence interval. Unit is one person-year. |
|
255 |
#' |
|
256 |
#' @examples |
|
257 |
#' h_incidence_rate_byar(200, 2) |
|
258 |
#' |
|
259 |
#' @export |
|
260 |
h_incidence_rate_byar <- function(person_years, |
|
261 |
n_events, |
|
262 |
alpha = 0.05) { |
|
263 | 1x |
checkmate::assert_number(person_years) |
264 | 1x |
checkmate::assert_number(n_events) |
265 | 1x |
assert_proportion_value(alpha) |
266 | ||
267 | 1x |
est <- n_events / person_years |
268 | 1x |
seg_1 <- n_events + 0.5 |
269 | 1x |
seg_2 <- 1 - 1 / (9 * (n_events + 0.5)) |
270 | 1x |
seg_3 <- stats::qnorm(1 - alpha / 2) * sqrt(1 / (n_events + 0.5)) / 3 |
271 | 1x |
lcl <- seg_1 * ((seg_2 - seg_3)^3) / person_years |
272 | 1x |
ucl <- seg_1 * ((seg_2 + seg_3) ^ 3) / person_years # styler: off |
273 | ||
274 | 1x |
list(rate = est, rate_ci = c(lcl, ucl)) |
275 |
} |
|
276 | ||
277 |
#' @describeIn h_incidence_rate Helper function to estimate the incidence rate and |
|
278 |
#' associated confidence interval. |
|
279 |
#' |
|
280 |
#' |
|
281 |
#' @keywords internal |
|
282 |
h_incidence_rate <- function(person_years, |
|
283 |
n_events, |
|
284 |
control = control_incidence_rate()) { |
|
285 | 4x |
alpha <- 1 - control$conf_level |
286 | 4x |
est <- switch(control$conf_type, |
287 | 4x |
normal = h_incidence_rate_normal(person_years, n_events, alpha), |
288 | 4x |
normal_log = h_incidence_rate_normal_log(person_years, n_events, alpha), |
289 | 4x |
exact = h_incidence_rate_exact(person_years, n_events, alpha), |
290 | 4x |
byar = h_incidence_rate_byar(person_years, n_events, alpha) |
291 |
) |
|
292 | ||
293 | 4x |
num_pt_year <- control$num_pt_year |
294 | 4x |
list( |
295 | 4x |
rate = est$rate * num_pt_year, |
296 | 4x |
rate_ci = est$rate_ci * num_pt_year |
297 |
) |
|
298 |
} |
1 |
#' Kaplan-Meier Plot |
|
2 |
#' |
|
3 |
#' @description `r lifecycle::badge("stable")` |
|
4 |
#' |
|
5 |
#' From a survival model, a graphic is rendered along with tabulated annotation |
|
6 |
#' including the number of patient at risk at given time and the median survival |
|
7 |
#' per group. |
|
8 |
#' |
|
9 |
#' @inheritParams grid::gTree |
|
10 |
#' @inheritParams argument_convention |
|
11 |
#' @param df (`data.frame`)\cr data set containing all analysis variables. |
|
12 |
#' @param variables (named `list`)\cr variable names. Details are: |
|
13 |
#' * `tte` (`numeric`)\cr variable indicating time-to-event duration values. |
|
14 |
#' * `is_event` (`logical`)\cr event variable. `TRUE` if event, `FALSE` if time to event is censored. |
|
15 |
#' * `arm` (`factor`)\cr the treatment group variable. |
|
16 |
#' * `strat` (`character` or `NULL`)\cr variable names indicating stratification factors. |
|
17 |
#' @param control_surv (`list`)\cr parameters for comparison details, specified by using |
|
18 |
#' the helper function [control_surv_timepoint()]. Some possible parameter options are: |
|
19 |
#' * `conf_level` (`proportion`)\cr confidence level of the interval for survival rate. |
|
20 |
#' * `conf_type` (`string`)\cr `"plain"` (default), `"log"`, `"log-log"` for confidence interval type, |
|
21 |
#' see more in [survival::survfit()]. Note that the option "none" is no longer supported. |
|
22 |
#' @param xticks (`numeric`, `number`, or `NULL`)\cr numeric vector of ticks or single number with spacing |
|
23 |
#' between ticks on the x axis. If `NULL` (default), [labeling::extended()] is used to determine |
|
24 |
#' an optimal tick position on the x axis. |
|
25 |
#' @param yval (`string`)\cr value of y-axis. Options are `Survival` (default) and `Failure` probability. |
|
26 |
#' @param censor_show (`flag`)\cr whether to show censored. |
|
27 |
#' @param xlab (`string`)\cr label of x-axis. |
|
28 |
#' @param ylab (`string`)\cr label of y-axis. |
|
29 |
#' @param ylim (`vector` of `numeric`)\cr vector of length 2 containing lower and upper limits for the y-axis. |
|
30 |
#' If `NULL` (default), the minimum and maximum y-values displayed are used as limits. |
|
31 |
#' @param title (`string`)\cr title for plot. |
|
32 |
#' @param footnotes (`string`)\cr footnotes for plot. |
|
33 |
#' @param col (`character`)\cr lines colors. Length of a vector should be equal |
|
34 |
#' to number of strata from [survival::survfit()]. |
|
35 |
#' @param lty (`numeric`)\cr line type. Length of a vector should be equal |
|
36 |
#' to number of strata from [survival::survfit()]. |
|
37 |
#' @param lwd (`numeric`)\cr line width. Length of a vector should be equal |
|
38 |
#' to number of strata from [survival::survfit()]. |
|
39 |
#' @param pch (`numeric`, `string`)\cr value or character of points symbol to indicate censored cases. |
|
40 |
#' @param size (`numeric`)\cr size of censored point, a class of `unit`. |
|
41 |
#' @param max_time (`numeric`)\cr maximum value to show on X axis. Only data values less than or up to |
|
42 |
#' this threshold value will be plotted (defaults to `NULL`). |
|
43 |
#' @param font_size (`number`)\cr font size to be used. |
|
44 |
#' @param ci_ribbon (`flag`)\cr draw the confidence interval around the Kaplan-Meier curve. |
|
45 |
#' @param ggtheme (`theme`)\cr a graphical theme as provided by `ggplot2` to control outlook of the Kaplan-Meier curve. |
|
46 |
#' @param annot_at_risk (`flag`)\cr compute and add the annotation table reporting the number of patient at risk |
|
47 |
#' matching the main grid of the Kaplan-Meier curve. |
|
48 |
#' @param annot_at_risk_title (`flag`)\cr whether the "Patients at Risk" title should be added above the `annot_at_risk` |
|
49 |
#' table. Has no effect if `annot_at_risk` is `FALSE`. Defaults to `TRUE`. |
|
50 |
#' @param annot_surv_med (`flag`)\cr compute and add the annotation table on the Kaplan-Meier curve estimating the |
|
51 |
#' median survival time per group. |
|
52 |
#' @param annot_coxph (`flag`)\cr add the annotation table from a [survival::coxph()] model. |
|
53 |
#' @param annot_stats (`string`)\cr statistics annotations to add to the plot. Options are |
|
54 |
#' `median` (median survival follow-up time) and `min` (minimum survival follow-up time). |
|
55 |
#' @param annot_stats_vlines (`flag`)\cr add vertical lines corresponding to each of the statistics |
|
56 |
#' specified by `annot_stats`. If `annot_stats` is `NULL` no lines will be added. |
|
57 |
#' @param control_coxph_pw (`list`)\cr parameters for comparison details, specified by using |
|
58 |
#' the helper function [control_coxph()]. Some possible parameter options are: |
|
59 |
#' * `pval_method` (`string`)\cr p-value method for testing hazard ratio = 1. |
|
60 |
#' Default method is `"log-rank"`, can also be set to `"wald"` or `"likelihood"`. |
|
61 |
#' * `ties` (`string`)\cr method for tie handling. Default is `"efron"`, |
|
62 |
#' can also be set to `"breslow"` or `"exact"`. See more in [survival::coxph()] |
|
63 |
#' * `conf_level` (`proportion`)\cr confidence level of the interval for HR. |
|
64 |
#' @param position_coxph (`numeric`)\cr x and y positions for plotting [survival::coxph()] model. |
|
65 |
#' @param position_surv_med (`numeric`)\cr x and y positions for plotting annotation table estimating median survival |
|
66 |
#' time per group. |
|
67 |
#' @param width_annots (named `list` of `unit`s)\cr a named list of widths for annotation tables with names `surv_med` |
|
68 |
#' (median survival time table) and `coxph` ([survival::coxph()] model table), where each value is the width |
|
69 |
#' (in units) to implement when printing the annotation table. |
|
70 |
#' |
|
71 |
#' @return A `grob` of class `gTree`. |
|
72 |
#' |
|
73 |
#' @examples |
|
74 |
#' \donttest{ |
|
75 |
#' library(dplyr) |
|
76 |
#' library(ggplot2) |
|
77 |
#' library(survival) |
|
78 |
#' library(grid) |
|
79 |
#' library(nestcolor) |
|
80 |
#' |
|
81 |
#' df <- tern_ex_adtte %>% |
|
82 |
#' filter(PARAMCD == "OS") %>% |
|
83 |
#' mutate(is_event = CNSR == 0) |
|
84 |
#' variables <- list(tte = "AVAL", is_event = "is_event", arm = "ARMCD") |
|
85 |
#' |
|
86 |
#' # 1. Example - basic option |
|
87 |
#' |
|
88 |
#' res <- g_km(df = df, variables = variables) |
|
89 |
#' res <- g_km(df = df, variables = variables, yval = "Failure") |
|
90 |
#' res <- g_km( |
|
91 |
#' df = df, |
|
92 |
#' variables = variables, |
|
93 |
#' control_surv = control_surv_timepoint(conf_level = 0.9), |
|
94 |
#' col = c("grey25", "grey50", "grey75"), |
|
95 |
#' annot_at_risk_title = FALSE |
|
96 |
#' ) |
|
97 |
#' res <- g_km(df = df, variables = variables, ggtheme = theme_minimal()) |
|
98 |
#' res <- g_km(df = df, variables = variables, ggtheme = theme_minimal(), lty = 1:3) |
|
99 |
#' res <- g_km(df = df, variables = variables, max = 2000) |
|
100 |
#' res <- g_km( |
|
101 |
#' df = df, |
|
102 |
#' variables = variables, |
|
103 |
#' annot_stats = c("min", "median"), |
|
104 |
#' annot_stats_vlines = TRUE |
|
105 |
#' ) |
|
106 |
#' |
|
107 |
#' # 2. Example - Arrange several KM curve on a single graph device |
|
108 |
#' |
|
109 |
#' # 2.1 Use case: A general graph on the top, a zoom on the bottom. |
|
110 |
#' grid.newpage() |
|
111 |
#' lyt <- grid.layout(nrow = 2, ncol = 1) %>% |
|
112 |
#' viewport(layout = .) %>% |
|
113 |
#' pushViewport() |
|
114 |
#' |
|
115 |
#' res <- g_km( |
|
116 |
#' df = df, variables = variables, newpage = FALSE, annot_surv_med = FALSE, |
|
117 |
#' vp = viewport(layout.pos.row = 1, layout.pos.col = 1) |
|
118 |
#' ) |
|
119 |
#' res <- g_km( |
|
120 |
#' df = df, variables = variables, max = 1000, newpage = FALSE, annot_surv_med = FALSE, |
|
121 |
#' ggtheme = theme_dark(), |
|
122 |
#' vp = viewport(layout.pos.row = 2, layout.pos.col = 1) |
|
123 |
#' ) |
|
124 |
#' |
|
125 |
#' # 2.1 Use case: No annotations on top, annotated graph on bottom |
|
126 |
#' grid.newpage() |
|
127 |
#' lyt <- grid.layout(nrow = 2, ncol = 1) %>% |
|
128 |
#' viewport(layout = .) %>% |
|
129 |
#' pushViewport() |
|
130 |
#' |
|
131 |
#' res <- g_km( |
|
132 |
#' df = df, variables = variables, newpage = FALSE, |
|
133 |
#' annot_surv_med = FALSE, annot_at_risk = FALSE, |
|
134 |
#' vp = viewport(layout.pos.row = 1, layout.pos.col = 1) |
|
135 |
#' ) |
|
136 |
#' res <- g_km( |
|
137 |
#' df = df, variables = variables, max = 2000, newpage = FALSE, annot_surv_med = FALSE, |
|
138 |
#' annot_at_risk = TRUE, |
|
139 |
#' ggtheme = theme_dark(), |
|
140 |
#' vp = viewport(layout.pos.row = 2, layout.pos.col = 1) |
|
141 |
#' ) |
|
142 |
#' |
|
143 |
#' # Add annotation from a pairwise coxph analysis |
|
144 |
#' g_km( |
|
145 |
#' df = df, variables = variables, |
|
146 |
#' annot_coxph = TRUE |
|
147 |
#' ) |
|
148 |
#' |
|
149 |
#' # Change widths/sizes of surv_med and coxph annotation tables. |
|
150 |
#' g_km( |
|
151 |
#' df = df, variables = c(variables, list(strat = "SEX")), |
|
152 |
#' annot_coxph = TRUE, |
|
153 |
#' width_annots = list(surv_med = grid::unit(2, "in"), coxph = grid::unit(3, "in")) |
|
154 |
#' ) |
|
155 |
#' |
|
156 |
#' g_km( |
|
157 |
#' df = df, variables = c(variables, list(strat = "SEX")), |
|
158 |
#' font_size = 15, |
|
159 |
#' annot_coxph = TRUE, |
|
160 |
#' control_coxph = control_coxph(pval_method = "wald", ties = "exact", conf_level = 0.99), |
|
161 |
#' position_coxph = c(0.5, 0.5) |
|
162 |
#' ) |
|
163 |
#' |
|
164 |
#' # Change position of the treatment group annotation table. |
|
165 |
#' g_km( |
|
166 |
#' df = df, variables = c(variables, list(strat = "SEX")), |
|
167 |
#' font_size = 15, |
|
168 |
#' annot_coxph = TRUE, |
|
169 |
#' control_coxph = control_coxph(pval_method = "wald", ties = "exact", conf_level = 0.99), |
|
170 |
#' position_surv_med = c(1, 0.7) |
|
171 |
#' ) |
|
172 |
#' } |
|
173 |
#' |
|
174 |
#' @export |
|
175 |
g_km <- function(df, |
|
176 |
variables, |
|
177 |
control_surv = control_surv_timepoint(), |
|
178 |
col = NULL, |
|
179 |
lty = NULL, |
|
180 |
lwd = .5, |
|
181 |
censor_show = TRUE, |
|
182 |
pch = 3, |
|
183 |
size = 2, |
|
184 |
max_time = NULL, |
|
185 |
xticks = NULL, |
|
186 |
xlab = "Days", |
|
187 |
yval = c("Survival", "Failure"), |
|
188 |
ylab = paste(yval, "Probability"), |
|
189 |
ylim = NULL, |
|
190 |
title = NULL, |
|
191 |
footnotes = NULL, |
|
192 |
draw = TRUE, |
|
193 |
newpage = TRUE, |
|
194 |
gp = NULL, |
|
195 |
vp = NULL, |
|
196 |
name = NULL, |
|
197 |
font_size = 12, |
|
198 |
ci_ribbon = FALSE, |
|
199 |
ggtheme = nestcolor::theme_nest(), |
|
200 |
annot_at_risk = TRUE, |
|
201 |
annot_at_risk_title = TRUE, |
|
202 |
annot_surv_med = TRUE, |
|
203 |
annot_coxph = FALSE, |
|
204 |
annot_stats = NULL, |
|
205 |
annot_stats_vlines = FALSE, |
|
206 |
control_coxph_pw = control_coxph(), |
|
207 |
position_coxph = c(-0.03, -0.02), |
|
208 |
position_surv_med = c(0.95, 0.9), |
|
209 |
width_annots = list(surv_med = grid::unit(0.3, "npc"), coxph = grid::unit(0.4, "npc"))) { |
|
210 | 8x |
checkmate::assert_list(variables) |
211 | 8x |
checkmate::assert_subset(c("tte", "arm", "is_event"), names(variables)) |
212 | 8x |
checkmate::assert_string(title, null.ok = TRUE) |
213 | 8x |
checkmate::assert_string(footnotes, null.ok = TRUE) |
214 | 8x |
checkmate::assert_character(col, null.ok = TRUE) |
215 | 8x |
checkmate::assert_subset(annot_stats, c("median", "min")) |
216 | 8x |
checkmate::assert_logical(annot_stats_vlines) |
217 | 8x |
checkmate::assert_true(all(sapply(width_annots, grid::is.unit))) |
218 | ||
219 | 8x |
tte <- variables$tte |
220 | 8x |
is_event <- variables$is_event |
221 | 8x |
arm <- variables$arm |
222 | ||
223 | 8x |
assert_valid_factor(df[[arm]]) |
224 | 8x |
assert_df_with_variables(df, list(tte = tte, is_event = is_event, arm = arm)) |
225 | 8x |
checkmate::assert_logical(df[[is_event]], min.len = 1, any.missing = FALSE) |
226 | 8x |
checkmate::assert_numeric(df[[tte]], min.len = 1, any.missing = FALSE) |
227 | ||
228 | 8x |
armval <- as.character(unique(df[[arm]])) |
229 | 8x |
if (annot_coxph && length(armval) < 2) { |
230 | ! |
stop(paste( |
231 | ! |
"When `annot_coxph` = TRUE, `df` must contain at least 2 levels of `variables$arm`", |
232 | ! |
"in order to calculate the hazard ratio." |
233 |
)) |
|
234 | 8x |
} else if (length(armval) > 1) { |
235 | 8x |
armval <- NULL |
236 |
} |
|
237 | 8x |
yval <- match.arg(yval) |
238 | 8x |
formula <- stats::as.formula(paste0("survival::Surv(", tte, ", ", is_event, ") ~ ", arm)) |
239 | 8x |
fit_km <- survival::survfit( |
240 | 8x |
formula = formula, |
241 | 8x |
data = df, |
242 | 8x |
conf.int = control_surv$conf_level, |
243 | 8x |
conf.type = control_surv$conf_type |
244 |
) |
|
245 | 8x |
data_plot <- h_data_plot( |
246 | 8x |
fit_km = fit_km, |
247 | 8x |
armval = armval, |
248 | 8x |
max_time = max_time |
249 |
) |
|
250 | ||
251 | 8x |
xticks <- h_xticks(data = data_plot, xticks = xticks, max_time = max_time) |
252 | 8x |
gg <- h_ggkm( |
253 | 8x |
data = data_plot, |
254 | 8x |
censor_show = censor_show, |
255 | 8x |
pch = pch, |
256 | 8x |
size = size, |
257 | 8x |
xticks = xticks, |
258 | 8x |
xlab = xlab, |
259 | 8x |
yval = yval, |
260 | 8x |
ylab = ylab, |
261 | 8x |
ylim = ylim, |
262 | 8x |
title = title, |
263 | 8x |
footnotes = footnotes, |
264 | 8x |
max_time = max_time, |
265 | 8x |
lwd = lwd, |
266 | 8x |
lty = lty, |
267 | 8x |
col = col, |
268 | 8x |
ggtheme = ggtheme, |
269 | 8x |
ci_ribbon = ci_ribbon |
270 |
) |
|
271 | ||
272 | 8x |
if (!is.null(annot_stats)) { |
273 | ! |
if ("median" %in% annot_stats) { |
274 | ! |
fit_km_all <- survival::survfit( |
275 | ! |
formula = stats::as.formula(paste0("survival::Surv(", tte, ", ", is_event, ") ~ ", 1)), |
276 | ! |
data = df, |
277 | ! |
conf.int = control_surv$conf_level, |
278 | ! |
conf.type = control_surv$conf_type |
279 |
) |
|
280 | ! |
gg <- gg + |
281 | ! |
geom_text( |
282 | ! |
size = 8 / ggplot2::.pt, col = 1, |
283 | ! |
x = stats::median(fit_km_all) + 0.065 * max(data_plot$time), |
284 | ! |
y = ifelse(yval == "Survival", 0.62, 0.38), |
285 | ! |
label = paste("Median F/U:\n", round(stats::median(fit_km_all), 1), tolower(df$AVALU[1])) |
286 |
) |
|
287 | ! |
if (annot_stats_vlines) { |
288 | ! |
gg <- gg + |
289 | ! |
geom_segment(aes(x = stats::median(fit_km_all), xend = stats::median(fit_km_all), y = -Inf, yend = Inf), |
290 | ! |
linetype = 2, col = "darkgray" |
291 |
) |
|
292 |
} |
|
293 |
} |
|
294 | ! |
if ("min" %in% annot_stats) { |
295 | ! |
min_fu <- min(df[[tte]]) |
296 | ! |
gg <- gg + |
297 | ! |
geom_text( |
298 | ! |
size = 8 / ggplot2::.pt, col = 1, |
299 | ! |
x = min_fu + max(data_plot$time) * ifelse(yval == "Survival", 0.05, 0.07), |
300 | ! |
y = ifelse(yval == "Survival", 1.0, 0.05), |
301 | ! |
label = paste("Min. F/U:\n", round(min_fu, 1), tolower(df$AVALU[1])) |
302 |
) |
|
303 | ! |
if (annot_stats_vlines) { |
304 | ! |
gg <- gg + |
305 | ! |
geom_segment(aes(x = min_fu, xend = min_fu, y = Inf, yend = -Inf), linetype = 2, col = "darkgray") |
306 |
} |
|
307 |
} |
|
308 | ! |
gg <- gg + ggplot2::guides(fill = ggplot2::guide_legend(override.aes = list(shape = NA, label = ""))) |
309 |
} |
|
310 | ||
311 | 8x |
g_el <- h_decompose_gg(gg) |
312 | ||
313 | 8x |
if (annot_at_risk) { |
314 |
# This is the content of the table that will be below the graph. |
|
315 | 6x |
annot_tbl <- summary(fit_km, time = xticks) |
316 | 6x |
annot_tbl <- if (is.null(fit_km$strata)) { |
317 | ! |
data.frame( |
318 | ! |
n.risk = annot_tbl$n.risk, |
319 | ! |
time = annot_tbl$time, |
320 | ! |
strata = as.factor(armval) |
321 |
) |
|
322 |
} else { |
|
323 | 6x |
strata_lst <- strsplit(sub("=", "equals", levels(annot_tbl$strata)), "equals") |
324 | 6x |
levels(annot_tbl$strata) <- matrix(unlist(strata_lst), ncol = 2, byrow = TRUE)[, 2] |
325 | 6x |
data.frame( |
326 | 6x |
n.risk = annot_tbl$n.risk, |
327 | 6x |
time = annot_tbl$time, |
328 | 6x |
strata = annot_tbl$strata |
329 |
) |
|
330 |
} |
|
331 | ||
332 | 6x |
grobs_patient <- h_grob_tbl_at_risk( |
333 | 6x |
data = data_plot, |
334 | 6x |
annot_tbl = annot_tbl, |
335 | 6x |
xlim = max(max_time, data_plot$time, xticks), |
336 | 6x |
title = annot_at_risk_title |
337 |
) |
|
338 |
} |
|
339 | ||
340 | 8x |
if (annot_at_risk || annot_surv_med || annot_coxph) { |
341 | 6x |
lyt <- h_km_layout( |
342 | 6x |
data = data_plot, g_el = g_el, title = title, footnotes = footnotes, |
343 | 6x |
annot_at_risk = annot_at_risk, annot_at_risk_title = annot_at_risk_title |
344 |
) |
|
345 | 6x |
at_risk_ttl <- as.numeric(annot_at_risk_title) |
346 | 6x |
ttl_row <- as.numeric(!is.null(title)) |
347 | 6x |
foot_row <- as.numeric(!is.null(footnotes)) |
348 | 6x |
km_grob <- grid::gTree( |
349 | 6x |
vp = grid::viewport(layout = lyt, height = .95, width = .95), |
350 | 6x |
children = grid::gList( |
351 |
# Title. |
|
352 | 6x |
if (ttl_row == 1) { |
353 | 1x |
grid::gTree( |
354 | 1x |
vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 2), |
355 | 1x |
children = grid::gList(grid::textGrob(label = title, x = grid::unit(0, "npc"), hjust = 0)) |
356 |
) |
|
357 |
}, |
|
358 | ||
359 |
# The Kaplan - Meier curve (top-right corner). |
|
360 | 6x |
grid::gTree( |
361 | 6x |
vp = grid::viewport(layout.pos.row = 1 + ttl_row, layout.pos.col = 2), |
362 | 6x |
children = grid::gList(g_el$panel) |
363 |
), |
|
364 | ||
365 |
# Survfit summary table (top-right corner). |
|
366 | 6x |
if (annot_surv_med) { |
367 | 5x |
grid::gTree( |
368 | 5x |
vp = grid::viewport(layout.pos.row = 1 + ttl_row, layout.pos.col = 2), |
369 | 5x |
children = h_grob_median_surv( |
370 | 5x |
fit_km = fit_km, |
371 | 5x |
armval = armval, |
372 | 5x |
x = position_surv_med[1], |
373 | 5x |
y = position_surv_med[2], |
374 | 5x |
width = if (!is.null(width_annots[["surv_med"]])) width_annots[["surv_med"]] else grid::unit(0.3, "npc"), |
375 | 5x |
ttheme = gridExtra::ttheme_default(base_size = font_size) |
376 |
) |
|
377 |
) |
|
378 |
}, |
|
379 | 6x |
if (annot_coxph) { |
380 | 1x |
grid::gTree( |
381 | 1x |
vp = grid::viewport(layout.pos.row = 1 + ttl_row, layout.pos.col = 2), |
382 | 1x |
children = h_grob_coxph( |
383 | 1x |
df = df, |
384 | 1x |
variables = variables, |
385 | 1x |
control_coxph_pw = control_coxph_pw, |
386 | 1x |
x = position_coxph[1], |
387 | 1x |
y = position_coxph[2], |
388 | 1x |
width = if (!is.null(width_annots[["coxph"]])) width_annots[["coxph"]] else grid::unit(0.4, "npc"), |
389 | 1x |
ttheme = gridExtra::ttheme_default( |
390 | 1x |
base_size = font_size, |
391 | 1x |
padding = grid::unit(c(1, .5), "lines"), |
392 | 1x |
core = list(bg_params = list(fill = c("grey95", "grey90"), alpha = .5)) |
393 |
) |
|
394 |
) |
|
395 |
) |
|
396 |
}, |
|
397 | ||
398 |
# Add the y-axis annotation (top-left corner). |
|
399 | 6x |
grid::gTree( |
400 | 6x |
vp = grid::viewport(layout.pos.row = 1 + ttl_row, layout.pos.col = 1), |
401 | 6x |
children = h_grob_y_annot(ylab = g_el$ylab, yaxis = g_el$yaxis) |
402 |
), |
|
403 | ||
404 |
# Add the x-axis annotation (second row below the Kaplan Meier Curve). |
|
405 | 6x |
grid::gTree( |
406 | 6x |
vp = grid::viewport(layout.pos.row = 2 + ttl_row, layout.pos.col = 2), |
407 | 6x |
children = grid::gList(rbind(g_el$xaxis, g_el$xlab)) |
408 |
), |
|
409 | ||
410 |
# Add the legend. |
|
411 | 6x |
grid::gTree( |
412 | 6x |
vp = grid::viewport(layout.pos.row = 3 + ttl_row, layout.pos.col = 2), |
413 | 6x |
children = grid::gList(g_el$guide) |
414 |
), |
|
415 | ||
416 |
# Add the table with patient-at-risk numbers. |
|
417 | 6x |
if (annot_at_risk && annot_at_risk_title) { |
418 | 6x |
grid::gTree( |
419 | 6x |
vp = grid::viewport(layout.pos.row = 4 + ttl_row, layout.pos.col = 1), |
420 | 6x |
children = grobs_patient$title |
421 |
) |
|
422 |
}, |
|
423 | 6x |
if (annot_at_risk) { |
424 | 6x |
grid::gTree( |
425 | 6x |
vp = grid::viewport(layout.pos.row = 4 + at_risk_ttl + ttl_row, layout.pos.col = 2), |
426 | 6x |
children = grobs_patient$at_risk |
427 |
) |
|
428 |
}, |
|
429 | 6x |
if (annot_at_risk) { |
430 | 6x |
grid::gTree( |
431 | 6x |
vp = grid::viewport(layout.pos.row = 4 + at_risk_ttl + ttl_row, layout.pos.col = 1), |
432 | 6x |
children = grobs_patient$label |
433 |
) |
|
434 |
}, |
|
435 | 6x |
if (annot_at_risk) { |
436 |
# Add the x-axis for the table. |
|
437 | 6x |
grid::gTree( |
438 | 6x |
vp = grid::viewport(layout.pos.row = 5 + at_risk_ttl + ttl_row, layout.pos.col = 2), |
439 | 6x |
children = grid::gList(rbind(g_el$xaxis, g_el$xlab)) |
440 |
) |
|
441 |
}, |
|
442 | ||
443 |
# Footnotes. |
|
444 | 6x |
if (foot_row == 1) { |
445 | 1x |
grid::gTree( |
446 | 1x |
vp = grid::viewport( |
447 | 1x |
layout.pos.row = ifelse(annot_at_risk, 6 + at_risk_ttl + ttl_row, 4 + ttl_row), |
448 | 1x |
layout.pos.col = 2 |
449 |
), |
|
450 | 1x |
children = grid::gList(grid::textGrob(label = footnotes, x = grid::unit(0, "npc"), hjust = 0)) |
451 |
) |
|
452 |
} |
|
453 |
) |
|
454 |
) |
|
455 | ||
456 | 6x |
result <- grid::gTree( |
457 | 6x |
vp = vp, |
458 | 6x |
gp = gp, |
459 | 6x |
name = name, |
460 | 6x |
children = grid::gList(km_grob) |
461 |
) |
|
462 |
} else { |
|
463 | 2x |
result <- grid::gTree( |
464 | 2x |
vp = vp, |
465 | 2x |
gp = gp, |
466 | 2x |
name = name, |
467 | 2x |
children = grid::gList(ggplot2::ggplotGrob(gg)) |
468 |
) |
|
469 |
} |
|
470 | ||
471 | 5x |
if (newpage && draw) grid::grid.newpage() |
472 | 5x |
if (draw) grid::grid.draw(result) |
473 | 8x |
invisible(result) |
474 |
} |
|
475 | ||
476 |
#' Helper function: tidy survival fit |
|
477 |
#' |
|
478 |
#' @description `r lifecycle::badge("stable")` |
|
479 |
#' |
|
480 |
#' Convert the survival fit data into a data frame designed for plotting |
|
481 |
#' within `g_km`. |
|
482 |
#' |
|
483 |
#' This starts from the [broom::tidy()] result, and then: |
|
484 |
#' * Post-processes the `strata` column into a factor. |
|
485 |
#' * Extends each stratum by an additional first row with time 0 and probability 1 so that |
|
486 |
#' downstream plot lines start at those coordinates. |
|
487 |
#' * Adds a `censor` column. |
|
488 |
#' * Filters the rows before `max_time`. |
|
489 |
#' |
|
490 |
#' @inheritParams g_km |
|
491 |
#' @param fit_km (`survfit`)\cr result of [survival::survfit()]. |
|
492 |
#' @param armval (`string`)\cr used as strata name when treatment arm variable only has one level. Default is `"All"`. |
|
493 |
#' |
|
494 |
#' @return A `tibble` with columns `time`, `n.risk`, `n.event`, `n.censor`, `estimate`, `std.error`, `conf.high`, |
|
495 |
#' `conf.low`, `strata`, and `censor`. |
|
496 |
#' |
|
497 |
#' @examples |
|
498 |
#' \donttest{ |
|
499 |
#' library(dplyr) |
|
500 |
#' library(survival) |
|
501 |
#' |
|
502 |
#' # Test with multiple arms |
|
503 |
#' tern_ex_adtte %>% |
|
504 |
#' filter(PARAMCD == "OS") %>% |
|
505 |
#' survfit(form = Surv(AVAL, 1 - CNSR) ~ ARMCD, data = .) %>% |
|
506 |
#' h_data_plot() |
|
507 |
#' |
|
508 |
#' # Test with single arm |
|
509 |
#' tern_ex_adtte %>% |
|
510 |
#' filter(PARAMCD == "OS", ARMCD == "ARM B") %>% |
|
511 |
#' survfit(form = Surv(AVAL, 1 - CNSR) ~ ARMCD, data = .) %>% |
|
512 |
#' h_data_plot(armval = "ARM B") |
|
513 |
#' } |
|
514 |
#' |
|
515 |
#' @export |
|
516 |
h_data_plot <- function(fit_km, |
|
517 |
armval = "All", |
|
518 |
max_time = NULL) { |
|
519 | 15x |
y <- broom::tidy(fit_km) |
520 | ||
521 | 15x |
if (!is.null(fit_km$strata)) { |
522 | 15x |
fit_km_var_level <- strsplit(sub("=", "equals", names(fit_km$strata)), "equals") |
523 | 15x |
strata_levels <- vapply(fit_km_var_level, FUN = "[", FUN.VALUE = "a", i = 2) |
524 | 15x |
strata_var_level <- strsplit(sub("=", "equals", y$strata), "equals") |
525 | 15x |
y$strata <- factor( |
526 | 15x |
vapply(strata_var_level, FUN = "[", FUN.VALUE = "a", i = 2), |
527 | 15x |
levels = strata_levels |
528 |
) |
|
529 |
} else { |
|
530 | ! |
y$strata <- armval |
531 |
} |
|
532 | ||
533 | 15x |
y_by_strata <- split(y, y$strata) |
534 | 15x |
y_by_strata_extended <- lapply( |
535 | 15x |
y_by_strata, |
536 | 15x |
FUN = function(tbl) { |
537 | 44x |
first_row <- tbl[1L, ] |
538 | 44x |
first_row$time <- 0 |
539 | 44x |
first_row$n.risk <- sum(first_row[, c("n.risk", "n.event", "n.censor")]) |
540 | 44x |
first_row$n.event <- first_row$n.censor <- 0 |
541 | 44x |
first_row$estimate <- first_row$conf.high <- first_row$conf.low <- 1 |
542 | 44x |
first_row$std.error <- 0 |
543 | 44x |
rbind( |
544 | 44x |
first_row, |
545 | 44x |
tbl |
546 |
) |
|
547 |
} |
|
548 |
) |
|
549 | 15x |
y <- do.call(rbind, y_by_strata_extended) |
550 | ||
551 | 15x |
y$censor <- ifelse(y$n.censor > 0, y$estimate, NA) |
552 | 15x |
if (!is.null(max_time)) { |
553 | 3x |
y <- y[y$time <= max(max_time), ] |
554 |
} |
|
555 | 15x |
y |
556 |
} |
|
557 | ||
558 |
#' Helper function: x tick positions |
|
559 |
#' |
|
560 |
#' @description `r lifecycle::badge("stable")` |
|
561 |
#' |
|
562 |
#' Calculate the positions of ticks on the x-axis. However, if `xticks` already |
|
563 |
#' exists it is kept as is. It is based on the same function `ggplot2` relies on, |
|
564 |
#' and is required in the graphic and the patient-at-risk annotation table. |
|
565 |
#' |
|
566 |
#' @inheritParams g_km |
|
567 |
#' @inheritParams h_ggkm |
|
568 |
#' |
|
569 |
#' @return A vector of positions to use for x-axis ticks on a `ggplot` object. |
|
570 |
#' |
|
571 |
#' @examples |
|
572 |
#' \donttest{ |
|
573 |
#' library(dplyr) |
|
574 |
#' library(survival) |
|
575 |
#' |
|
576 |
#' data <- tern_ex_adtte %>% |
|
577 |
#' filter(PARAMCD == "OS") %>% |
|
578 |
#' survfit(form = Surv(AVAL, 1 - CNSR) ~ ARMCD, data = .) %>% |
|
579 |
#' h_data_plot() |
|
580 |
#' |
|
581 |
#' h_xticks(data) |
|
582 |
#' h_xticks(data, xticks = seq(0, 3000, 500)) |
|
583 |
#' h_xticks(data, xticks = 500) |
|
584 |
#' h_xticks(data, xticks = 500, max_time = 6000) |
|
585 |
#' h_xticks(data, xticks = c(0, 500), max_time = 300) |
|
586 |
#' h_xticks(data, xticks = 500, max_time = 300) |
|
587 |
#' } |
|
588 |
#' |
|
589 |
#' @export |
|
590 |
h_xticks <- function(data, xticks = NULL, max_time = NULL) { |
|
591 | 15x |
if (is.null(xticks)) { |
592 | 9x |
if (is.null(max_time)) { |
593 | 7x |
labeling::extended(range(data$time)[1], range(data$time)[2], m = 5) |
594 |
} else { |
|
595 | 2x |
labeling::extended(range(data$time)[1], max(range(data$time)[2], max_time), m = 5) |
596 |
} |
|
597 | 6x |
} else if (checkmate::test_number(xticks)) { |
598 | 3x |
if (is.null(max_time)) { |
599 | 2x |
seq(0, max(data$time), xticks) |
600 |
} else { |
|
601 | 1x |
seq(0, max(data$time, max_time), xticks) |
602 |
} |
|
603 | 3x |
} else if (is.numeric(xticks)) { |
604 | 2x |
xticks |
605 |
} else { |
|
606 | 1x |
stop( |
607 | 1x |
paste( |
608 | 1x |
"xticks should be either `NULL`", |
609 | 1x |
"or a single number (interval between x ticks)", |
610 | 1x |
"or a numeric vector (position of ticks on the x axis)" |
611 |
) |
|
612 |
) |
|
613 |
} |
|
614 |
} |
|
615 | ||
616 |
#' Helper function: KM plot |
|
617 |
#' |
|
618 |
#' @description `r lifecycle::badge("stable")` |
|
619 |
#' |
|
620 |
#' Draw the Kaplan-Meier plot using `ggplot2`. |
|
621 |
#' |
|
622 |
#' @inheritParams g_km |
|
623 |
#' @param data (`data.frame`)\cr survival data as pre-processed by `h_data_plot`. |
|
624 |
#' |
|
625 |
#' @return A `ggplot` object. |
|
626 |
#' |
|
627 |
#' @examples |
|
628 |
#' \donttest{ |
|
629 |
#' library(dplyr) |
|
630 |
#' library(survival) |
|
631 |
#' |
|
632 |
#' fit_km <- tern_ex_adtte %>% |
|
633 |
#' filter(PARAMCD == "OS") %>% |
|
634 |
#' survfit(form = Surv(AVAL, 1 - CNSR) ~ ARMCD, data = .) |
|
635 |
#' data_plot <- h_data_plot(fit_km = fit_km) |
|
636 |
#' xticks <- h_xticks(data = data_plot) |
|
637 |
#' gg <- h_ggkm( |
|
638 |
#' data = data_plot, |
|
639 |
#' censor_show = TRUE, |
|
640 |
#' xticks = xticks, |
|
641 |
#' xlab = "Days", |
|
642 |
#' yval = "Survival", |
|
643 |
#' ylab = "Survival Probability", |
|
644 |
#' title = "Survival" |
|
645 |
#' ) |
|
646 |
#' gg |
|
647 |
#' } |
|
648 |
#' |
|
649 |
#' @export |
|
650 |
h_ggkm <- function(data, |
|
651 |
xticks = NULL, |
|
652 |
yval = "Survival", |
|
653 |
censor_show, |
|
654 |
xlab, |
|
655 |
ylab, |
|
656 |
ylim = NULL, |
|
657 |
title, |
|
658 |
footnotes = NULL, |
|
659 |
max_time = NULL, |
|
660 |
lwd = 1, |
|
661 |
lty = NULL, |
|
662 |
pch = 3, |
|
663 |
size = 2, |
|
664 |
col = NULL, |
|
665 |
ci_ribbon = FALSE, |
|
666 |
ggtheme = nestcolor::theme_nest()) { |
|
667 | 8x |
checkmate::assert_numeric(lty, null.ok = TRUE) |
668 | 8x |
checkmate::assert_character(col, null.ok = TRUE) |
669 | ||
670 | 8x |
if (is.null(ylim)) { |
671 | 8x |
data_lims <- data |
672 | 1x |
if (yval == "Failure") data_lims[["estimate"]] <- 1 - data_lims[["estimate"]] |
673 | 8x |
if (!is.null(max_time)) { |
674 | 1x |
y_lwr <- min(data_lims[data_lims$time < max_time, ][["estimate"]]) |
675 | 1x |
y_upr <- max(data_lims[data_lims$time < max_time, ][["estimate"]]) |
676 |
} else { |
|
677 | 7x |
y_lwr <- min(data_lims[["estimate"]]) |
678 | 7x |
y_upr <- max(data_lims[["estimate"]]) |
679 |
} |
|
680 | 8x |
ylim <- c(y_lwr, y_upr) |
681 |
} |
|
682 | 8x |
checkmate::assert_numeric(ylim, finite = TRUE, any.missing = FALSE, len = 2, sorted = TRUE) |
683 | ||
684 |
# change estimates of survival to estimates of failure (1 - survival) |
|
685 | 8x |
if (yval == "Failure") { |
686 | 1x |
data$estimate <- 1 - data$estimate |
687 | 1x |
data[c("conf.high", "conf.low")] <- list(1 - data$conf.low, 1 - data$conf.high) |
688 | 1x |
data$censor <- 1 - data$censor |
689 |
} |
|
690 | ||
691 | 8x |
gg <- { |
692 | 8x |
ggplot2::ggplot( |
693 | 8x |
data = data, |
694 | 8x |
mapping = ggplot2::aes( |
695 | 8x |
x = .data[["time"]], |
696 | 8x |
y = .data[["estimate"]], |
697 | 8x |
ymin = .data[["conf.low"]], |
698 | 8x |
ymax = .data[["conf.high"]], |
699 | 8x |
color = .data[["strata"]], |
700 | 8x |
fill = .data[["strata"]] |
701 |
) |
|
702 |
) + |
|
703 | 8x |
ggplot2::geom_hline(yintercept = 0) |
704 |
} |
|
705 | ||
706 | 8x |
if (ci_ribbon) { |
707 | 1x |
gg <- gg + ggplot2::geom_ribbon(alpha = .3, lty = 0) |
708 |
} |
|
709 | ||
710 | 8x |
gg <- if (is.null(lty)) { |
711 | 7x |
gg + |
712 | 7x |
ggplot2::geom_step(linewidth = lwd) |
713 | 8x |
} else if (checkmate::test_number(lty)) { |
714 | 1x |
gg + |
715 | 1x |
ggplot2::geom_step(linewidth = lwd, lty = lty) |
716 | 8x |
} else if (is.numeric(lty)) { |
717 | ! |
gg + |
718 | ! |
ggplot2::geom_step(mapping = ggplot2::aes(linetype = .data[["strata"]]), linewidth = lwd) + |
719 | ! |
ggplot2::scale_linetype_manual(values = lty) |
720 |
} |
|
721 | ||
722 | 8x |
gg <- gg + |
723 | 8x |
ggplot2::coord_cartesian(ylim = ylim) + |
724 | 8x |
ggplot2::labs(x = xlab, y = ylab, title = title, caption = footnotes) |
725 | ||
726 | 8x |
if (!is.null(col)) { |
727 | ! |
gg <- gg + |
728 | ! |
ggplot2::scale_color_manual(values = col) + |
729 | ! |
ggplot2::scale_fill_manual(values = col) |
730 |
} |
|
731 | 8x |
if (censor_show) { |
732 | 8x |
dt <- data[data$n.censor != 0, ] |
733 | 8x |
dt$censor_lbl <- factor("Censored") |
734 | ||
735 | 8x |
gg <- gg + ggplot2::geom_point( |
736 | 8x |
data = dt, |
737 | 8x |
ggplot2::aes( |
738 | 8x |
x = .data[["time"]], |
739 | 8x |
y = .data[["censor"]], |
740 | 8x |
shape = .data[["censor_lbl"]] |
741 |
), |
|
742 | 8x |
size = size, |
743 | 8x |
show.legend = TRUE, |
744 | 8x |
inherit.aes = TRUE |
745 |
) + |
|
746 | 8x |
ggplot2::scale_shape_manual(name = NULL, values = pch) + |
747 | 8x |
ggplot2::guides( |
748 | 8x |
shape = ggplot2::guide_legend(override.aes = list(linetype = NA)), |
749 | 8x |
fill = ggplot2::guide_legend(override.aes = list(shape = NA)) |
750 |
) |
|
751 |
} |
|
752 | ||
753 | 8x |
if (!is.null(max_time) && !is.null(xticks)) { |
754 | 1x |
gg <- gg + ggplot2::scale_x_continuous(breaks = xticks, limits = c(min(0, xticks), max(c(xticks, max_time)))) |
755 | 7x |
} else if (!is.null(xticks)) { |
756 | 7x |
if (max(data$time) <= max(xticks)) { |
757 | 6x |
gg <- gg + ggplot2::scale_x_continuous(breaks = xticks, limits = c(min(0, min(xticks)), max(xticks))) |
758 |
} else { |
|
759 | 1x |
gg <- gg + ggplot2::scale_x_continuous(breaks = xticks) |
760 |
} |
|
761 | ! |
} else if (!is.null(max_time)) { |
762 | ! |
gg <- gg + ggplot2::scale_x_continuous(limits = c(0, max_time)) |
763 |
} |
|
764 | ||
765 | 8x |
if (!is.null(ggtheme)) { |
766 | 8x |
gg <- gg + ggtheme |
767 |
} |
|
768 | ||
769 | 8x |
gg + ggplot2::theme( |
770 | 8x |
legend.position = "bottom", |
771 | 8x |
legend.title = ggplot2::element_blank(), |
772 | 8x |
legend.key.height = unit(0.02, "npc"), |
773 | 8x |
panel.grid.major.x = ggplot2::element_line(linewidth = 2) |
774 |
) |
|
775 |
} |
|
776 | ||
777 |
#' `ggplot` Decomposition |
|
778 |
#' |
|
779 |
#' @description `r lifecycle::badge("stable")` |
|
780 |
#' |
|
781 |
#' The elements composing the `ggplot` are extracted and organized in a `list`. |
|
782 |
#' |
|
783 |
#' @param gg (`ggplot`)\cr a graphic to decompose. |
|
784 |
#' |
|
785 |
#' @return A named `list` with elements: |
|
786 |
#' * `panel`: The panel. |
|
787 |
#' * `yaxis`: The y-axis. |
|
788 |
#' * `xaxis`: The x-axis. |
|
789 |
#' * `xlab`: The x-axis label. |
|
790 |
#' * `ylab`: The y-axis label. |
|
791 |
#' * `guide`: The legend. |
|
792 |
#' |
|
793 |
#' @examples |
|
794 |
#' \donttest{ |
|
795 |
#' library(dplyr) |
|
796 |
#' library(survival) |
|
797 |
#' library(grid) |
|
798 |
#' |
|
799 |
#' fit_km <- tern_ex_adtte %>% |
|
800 |
#' filter(PARAMCD == "OS") %>% |
|
801 |
#' survfit(form = Surv(AVAL, 1 - CNSR) ~ ARMCD, data = .) |
|
802 |
#' data_plot <- h_data_plot(fit_km = fit_km) |
|
803 |
#' xticks <- h_xticks(data = data_plot) |
|
804 |
#' gg <- h_ggkm( |
|
805 |
#' data = data_plot, |
|
806 |
#' yval = "Survival", |
|
807 |
#' censor_show = TRUE, |
|
808 |
#' xticks = xticks, xlab = "Days", ylab = "Survival Probability", |
|
809 |
#' title = "tt", |
|
810 |
#' footnotes = "ff" |
|
811 |
#' ) |
|
812 |
#' |
|
813 |
#' g_el <- h_decompose_gg(gg) |
|
814 |
#' grid::grid.newpage() |
|
815 |
#' grid.rect(gp = grid::gpar(lty = 1, col = "red", fill = "gray85", lwd = 5)) |
|
816 |
#' grid::grid.draw(g_el$panel) |
|
817 |
#' |
|
818 |
#' grid::grid.newpage() |
|
819 |
#' grid.rect(gp = grid::gpar(lty = 1, col = "royalblue", fill = "gray85", lwd = 5)) |
|
820 |
#' grid::grid.draw(with(g_el, cbind(ylab, yaxis))) |
|
821 |
#' } |
|
822 |
#' |
|
823 |
#' @export |
|
824 |
h_decompose_gg <- function(gg) { |
|
825 | 8x |
g_el <- ggplot2::ggplotGrob(gg) |
826 | 8x |
y <- c( |
827 | 8x |
panel = "panel", |
828 | 8x |
yaxis = "axis-l", |
829 | 8x |
xaxis = "axis-b", |
830 | 8x |
xlab = "xlab-b", |
831 | 8x |
ylab = "ylab-l", |
832 | 8x |
guide = "guide" |
833 |
) |
|
834 | 8x |
lapply(X = y, function(x) gtable::gtable_filter(g_el, x)) |
835 |
} |
|
836 | ||
837 |
#' Helper: KM Layout |
|
838 |
#' |
|
839 |
#' @description `r lifecycle::badge("stable")` |
|
840 |
#' |
|
841 |
#' Prepares a (5 rows) x (2 cols) layout for the Kaplan-Meier curve. |
|
842 |
#' |
|
843 |
#' @inheritParams g_km |
|
844 |
#' @inheritParams h_ggkm |
|
845 |
#' @param g_el (`list` of `gtable`)\cr list as obtained by `h_decompose_gg()`. |
|
846 |
#' @param annot_at_risk (`flag`)\cr compute and add the annotation table reporting the number of |
|
847 |
#' patient at risk matching the main grid of the Kaplan-Meier curve. |
|
848 |
#' |
|
849 |
#' @return A grid layout. |
|
850 |
#' |
|
851 |
#' @details The layout corresponds to a grid of two columns and five rows of unequal dimensions. Most of the |
|
852 |
#' dimension are fixed, only the curve is flexible and will accommodate with the remaining free space. |
|
853 |
#' * The left column gets the annotation of the `ggplot` (y-axis) and the names of the strata for the patient |
|
854 |
#' at risk tabulation. The main constraint is about the width of the columns which must allow the writing of |
|
855 |
#' the strata name. |
|
856 |
#' * The right column receive the `ggplot`, the legend, the x-axis and the patient at risk table. |
|
857 |
#' |
|
858 |
#' @examples |
|
859 |
#' \donttest{ |
|
860 |
#' library(dplyr) |
|
861 |
#' library(survival) |
|
862 |
#' library(grid) |
|
863 |
#' |
|
864 |
#' fit_km <- tern_ex_adtte %>% |
|
865 |
#' filter(PARAMCD == "OS") %>% |
|
866 |
#' survfit(form = Surv(AVAL, 1 - CNSR) ~ ARMCD, data = .) |
|
867 |
#' data_plot <- h_data_plot(fit_km = fit_km) |
|
868 |
#' xticks <- h_xticks(data = data_plot) |
|
869 |
#' gg <- h_ggkm( |
|
870 |
#' data = data_plot, |
|
871 |
#' censor_show = TRUE, |
|
872 |
#' xticks = xticks, xlab = "Days", ylab = "Survival Probability", |
|
873 |
#' title = "tt", footnotes = "ff", yval = "Survival" |
|
874 |
#' ) |
|
875 |
#' g_el <- h_decompose_gg(gg) |
|
876 |
#' lyt <- h_km_layout(data = data_plot, g_el = g_el, title = "t", footnotes = "f") |
|
877 |
#' grid.show.layout(lyt) |
|
878 |
#' } |
|
879 |
#' |
|
880 |
#' @export |
|
881 |
h_km_layout <- function(data, g_el, title, footnotes, annot_at_risk = TRUE, annot_at_risk_title = TRUE) { |
|
882 | 6x |
txtlines <- levels(as.factor(data$strata)) |
883 | 6x |
nlines <- nlevels(as.factor(data$strata)) |
884 | 6x |
col_annot_width <- max( |
885 | 6x |
c( |
886 | 6x |
as.numeric(grid::convertX(g_el$yaxis$width + g_el$ylab$width, "pt")), |
887 | 6x |
as.numeric( |
888 | 6x |
grid::convertX( |
889 | 6x |
grid::stringWidth(txtlines) + grid::unit(7, "pt"), "pt" |
890 |
) |
|
891 |
) |
|
892 |
) |
|
893 |
) |
|
894 | ||
895 | 6x |
ttl_row <- as.numeric(!is.null(title)) |
896 | 6x |
foot_row <- as.numeric(!is.null(footnotes)) |
897 | 6x |
no_tbl_ind <- c() |
898 | 6x |
ht_x <- c() |
899 | 6x |
ht_units <- c() |
900 | ||
901 | 6x |
if (ttl_row == 1) { |
902 | 1x |
no_tbl_ind <- c(no_tbl_ind, TRUE) |
903 | 1x |
ht_x <- c(ht_x, 2) |
904 | 1x |
ht_units <- c(ht_units, "lines") |
905 |
} |
|
906 | ||
907 | 6x |
no_tbl_ind <- c(no_tbl_ind, rep(TRUE, 3), rep(FALSE, 2)) |
908 | 6x |
ht_x <- c( |
909 | 6x |
ht_x, |
910 | 6x |
1, |
911 | 6x |
grid::convertX(with(g_el, xaxis$height + ylab$width), "pt") + grid::unit(5, "pt"), |
912 | 6x |
grid::convertX(g_el$guide$heights, "pt") + grid::unit(2, "pt"), |
913 | 6x |
1, |
914 | 6x |
nlines + 0.5, |
915 | 6x |
grid::convertX(with(g_el, xaxis$height + ylab$width), "pt") |
916 |
) |
|
917 | 6x |
ht_units <- c( |
918 | 6x |
ht_units, |
919 | 6x |
"null", |
920 | 6x |
"pt", |
921 | 6x |
"pt", |
922 | 6x |
"lines", |
923 | 6x |
"lines", |
924 | 6x |
"pt" |
925 |
) |
|
926 | ||
927 | 6x |
if (foot_row == 1) { |
928 | 1x |
no_tbl_ind <- c(no_tbl_ind, TRUE) |
929 | 1x |
ht_x <- c(ht_x, 1) |
930 | 1x |
ht_units <- c(ht_units, "lines") |
931 |
} |
|
932 | 6x |
if (annot_at_risk) { |
933 | 6x |
no_at_risk_tbl <- rep(TRUE, 6 + ttl_row + foot_row) |
934 | 6x |
if (!annot_at_risk_title) { |
935 | ! |
no_at_risk_tbl[length(no_at_risk_tbl) - 2 - foot_row] <- FALSE |
936 |
} |
|
937 |
} else { |
|
938 | ! |
no_at_risk_tbl <- no_tbl_ind |
939 |
} |
|
940 | ||
941 | 6x |
grid::grid.layout( |
942 | 6x |
nrow = sum(no_at_risk_tbl), ncol = 2, |
943 | 6x |
widths = grid::unit(c(col_annot_width, 1), c("pt", "null")), |
944 | 6x |
heights = grid::unit( |
945 | 6x |
x = ht_x[no_at_risk_tbl], |
946 | 6x |
units = ht_units[no_at_risk_tbl] |
947 |
) |
|
948 |
) |
|
949 |
} |
|
950 | ||
951 |
#' Helper: Patient-at-Risk Grobs |
|
952 |
#' |
|
953 |
#' @description `r lifecycle::badge("stable")` |
|
954 |
#' |
|
955 |
#' Two graphical objects are obtained, one corresponding to row labeling and the second to the table of |
|
956 |
#' numbers of patients at risk. If `title = TRUE`, a third object corresponding to the table title is |
|
957 |
#' also obtained. |
|
958 |
#' |
|
959 |
#' @inheritParams g_km |
|
960 |
#' @inheritParams h_ggkm |
|
961 |
#' @param annot_tbl (`data.frame`)\cr annotation as prepared by [survival::summary.survfit()] which |
|
962 |
#' includes the number of patients at risk at given time points. |
|
963 |
#' @param xlim (`numeric`)\cr the maximum value on the x-axis (used to |
|
964 |
#' ensure the at risk table aligns with the KM graph). |
|
965 |
#' @param title (`flag`)\cr whether the "Patients at Risk" title should be added above the `annot_at_risk` |
|
966 |
#' table. Has no effect if `annot_at_risk` is `FALSE`. Defaults to `TRUE`. |
|
967 |
#' |
|
968 |
#' @return A named `list` of two `gTree` objects if `title = FALSE`: `at_risk` and `label`, or three |
|
969 |
#' `gTree` objects if `title = TRUE`: `at_risk`, `label`, and `title`. |
|
970 |
#' |
|
971 |
#' @examples |
|
972 |
#' \donttest{ |
|
973 |
#' library(dplyr) |
|
974 |
#' library(survival) |
|
975 |
#' library(grid) |
|
976 |
#' |
|
977 |
#' fit_km <- tern_ex_adtte %>% |
|
978 |
#' filter(PARAMCD == "OS") %>% |
|
979 |
#' survfit(form = Surv(AVAL, 1 - CNSR) ~ ARMCD, data = .) |
|
980 |
#' |
|
981 |
#' data_plot <- h_data_plot(fit_km = fit_km) |
|
982 |
#' |
|
983 |
#' xticks <- h_xticks(data = data_plot) |
|
984 |
#' |
|
985 |
#' gg <- h_ggkm( |
|
986 |
#' data = data_plot, |
|
987 |
#' censor_show = TRUE, |
|
988 |
#' xticks = xticks, xlab = "Days", ylab = "Survival Probability", |
|
989 |
#' title = "tt", footnotes = "ff", yval = "Survival" |
|
990 |
#' ) |
|
991 |
#' |
|
992 |
#' # The annotation table reports the patient at risk for a given strata and |
|
993 |
#' # time (`xticks`). |
|
994 |
#' annot_tbl <- summary(fit_km, time = xticks) |
|
995 |
#' if (is.null(fit_km$strata)) { |
|
996 |
#' annot_tbl <- with(annot_tbl, data.frame(n.risk = n.risk, time = time, strata = "All")) |
|
997 |
#' } else { |
|
998 |
#' strata_lst <- strsplit(sub("=", "equals", levels(annot_tbl$strata)), "equals") |
|
999 |
#' levels(annot_tbl$strata) <- matrix(unlist(strata_lst), ncol = 2, byrow = TRUE)[, 2] |
|
1000 |
#' annot_tbl <- data.frame( |
|
1001 |
#' n.risk = annot_tbl$n.risk, |
|
1002 |
#' time = annot_tbl$time, |
|
1003 |
#' strata = annot_tbl$strata |
|
1004 |
#' ) |
|
1005 |
#' } |
|
1006 |
#' |
|
1007 |
#' # The annotation table is transformed into a grob. |
|
1008 |
#' tbl <- h_grob_tbl_at_risk(data = data_plot, annot_tbl = annot_tbl, xlim = max(xticks)) |
|
1009 |
#' |
|
1010 |
#' # For the representation, the layout is estimated for which the decomposition |
|
1011 |
#' # of the graphic element is necessary. |
|
1012 |
#' g_el <- h_decompose_gg(gg) |
|
1013 |
#' lyt <- h_km_layout(data = data_plot, g_el = g_el, title = "t", footnotes = "f") |
|
1014 |
#' |
|
1015 |
#' grid::grid.newpage() |
|
1016 |
#' pushViewport(viewport(layout = lyt, height = .95, width = .95)) |
|
1017 |
#' grid.rect(gp = grid::gpar(lty = 1, col = "purple", fill = "gray85", lwd = 1)) |
|
1018 |
#' pushViewport(viewport(layout.pos.row = 3:4, layout.pos.col = 2)) |
|
1019 |
#' grid.rect(gp = grid::gpar(lty = 1, col = "orange", fill = "gray85", lwd = 1)) |
|
1020 |
#' grid::grid.draw(tbl$at_risk) |
|
1021 |
#' popViewport() |
|
1022 |
#' pushViewport(viewport(layout.pos.row = 3:4, layout.pos.col = 1)) |
|
1023 |
#' grid.rect(gp = grid::gpar(lty = 1, col = "green3", fill = "gray85", lwd = 1)) |
|
1024 |
#' grid::grid.draw(tbl$label) |
|
1025 |
#' } |
|
1026 |
#' |
|
1027 |
#' @export |
|
1028 |
h_grob_tbl_at_risk <- function(data, annot_tbl, xlim, title = TRUE) { |
|
1029 | 6x |
txtlines <- levels(as.factor(data$strata)) |
1030 | 6x |
nlines <- nlevels(as.factor(data$strata)) |
1031 | 6x |
y_int <- annot_tbl$time[2] - annot_tbl$time[1] |
1032 | 6x |
annot_tbl <- expand.grid( |
1033 | 6x |
time = seq(0, xlim, y_int), |
1034 | 6x |
strata = unique(annot_tbl$strata) |
1035 | 6x |
) %>% dplyr::left_join(annot_tbl, by = c("time", "strata")) |
1036 | 6x |
annot_tbl[is.na(annot_tbl)] <- 0 |
1037 | 6x |
y_str_unit <- as.numeric(annot_tbl$strata) |
1038 | 6x |
vp_table <- grid::plotViewport(margins = grid::unit(c(0, 0, 0, 0), "lines")) |
1039 | 6x |
if (title) { |
1040 | 6x |
gb_table_title <- grid::gList( |
1041 | 6x |
grid::textGrob( |
1042 | 6x |
label = "Patients at Risk:", |
1043 | 6x |
x = 1, |
1044 | 6x |
y = grid::unit(0.2, "native"), |
1045 | 6x |
gp = grid::gpar(fontface = "bold", fontsize = 10) |
1046 |
) |
|
1047 |
) |
|
1048 |
} |
|
1049 | 6x |
gb_table_left_annot <- grid::gList( |
1050 | 6x |
grid::rectGrob( |
1051 | 6x |
x = 0, y = grid::unit(c(1:nlines) - 1, "lines"), |
1052 | 6x |
gp = grid::gpar(fill = c("gray95", "gray90"), alpha = 1, col = "white"), |
1053 | 6x |
height = grid::unit(1, "lines"), just = "bottom", hjust = 0 |
1054 |
), |
|
1055 | 6x |
grid::textGrob( |
1056 | 6x |
label = unique(annot_tbl$strata), |
1057 | 6x |
x = 0.5, |
1058 | 6x |
y = grid::unit( |
1059 | 6x |
(max(unique(y_str_unit)) - unique(y_str_unit)) + 0.75, |
1060 | 6x |
"native" |
1061 |
), |
|
1062 | 6x |
gp = grid::gpar(fontface = "italic", fontsize = 10) |
1063 |
) |
|
1064 |
) |
|
1065 | 6x |
gb_patient_at_risk <- grid::gList( |
1066 | 6x |
grid::rectGrob( |
1067 | 6x |
x = 0, y = grid::unit(c(1:nlines) - 1, "lines"), |
1068 | 6x |
gp = grid::gpar(fill = c("gray95", "gray90"), alpha = 1, col = "white"), |
1069 | 6x |
height = grid::unit(1, "lines"), just = "bottom", hjust = 0 |
1070 |
), |
|
1071 | 6x |
grid::textGrob( |
1072 | 6x |
label = annot_tbl$n.risk, |
1073 | 6x |
x = grid::unit(annot_tbl$time, "native"), |
1074 | 6x |
y = grid::unit( |
1075 | 6x |
(max(y_str_unit) - y_str_unit) + .5, |
1076 | 6x |
"line" |
1077 | 6x |
) # maybe native |
1078 |
) |
|
1079 |
) |
|
1080 | ||
1081 | 6x |
ret <- list( |
1082 | 6x |
at_risk = grid::gList( |
1083 | 6x |
grid::gTree( |
1084 | 6x |
vp = vp_table, |
1085 | 6x |
children = grid::gList( |
1086 | 6x |
grid::gTree( |
1087 | 6x |
vp = grid::dataViewport( |
1088 | 6x |
xscale = c(0, xlim) + c(-0.05, 0.05) * xlim, |
1089 | 6x |
yscale = c(0, nlines + 1), |
1090 | 6x |
extension = c(0.05, 0) |
1091 |
), |
|
1092 | 6x |
children = grid::gList(gb_patient_at_risk) |
1093 |
) |
|
1094 |
) |
|
1095 |
) |
|
1096 |
), |
|
1097 | 6x |
label = grid::gList( |
1098 | 6x |
grid::gTree( |
1099 | 6x |
vp = grid::viewport(width = max(grid::stringWidth(txtlines))), |
1100 | 6x |
children = grid::gList( |
1101 | 6x |
grid::gTree( |
1102 | 6x |
vp = grid::dataViewport( |
1103 | 6x |
xscale = 0:1, |
1104 | 6x |
yscale = c(0, nlines + 1), |
1105 | 6x |
extension = c(0.0, 0) |
1106 |
), |
|
1107 | 6x |
children = grid::gList(gb_table_left_annot) |
1108 |
) |
|
1109 |
) |
|
1110 |
) |
|
1111 |
) |
|
1112 |
) |
|
1113 | ||
1114 | 6x |
if (title) { |
1115 | 6x |
ret[["title"]] <- grid::gList( |
1116 | 6x |
grid::gTree( |
1117 | 6x |
vp = grid::viewport(width = max(grid::stringWidth(txtlines))), |
1118 | 6x |
children = grid::gList( |
1119 | 6x |
grid::gTree( |
1120 | 6x |
vp = grid::dataViewport( |
1121 | 6x |
xscale = 0:1, |
1122 | 6x |
yscale = c(0, 1), |
1123 | 6x |
extension = c(0, 0) |
1124 |
), |
|
1125 | 6x |
children = grid::gList(gb_table_title) |
1126 |
) |
|
1127 |
) |
|
1128 |
) |
|
1129 |
) |
|
1130 |
} |
|
1131 | ||
1132 | 6x |
ret |
1133 |
} |
|
1134 | ||
1135 |
#' Helper Function: Survival Estimations |
|
1136 |
#' |
|
1137 |
#' @description `r lifecycle::badge("stable")` |
|
1138 |
#' |
|
1139 |
#' Transform a survival fit to a table with groups in rows characterized by N, median and confidence interval. |
|
1140 |
#' |
|
1141 |
#' @inheritParams h_data_plot |
|
1142 |
#' |
|
1143 |
#' @return A summary table with statistics `N`, `Median`, and `XX% CI` (`XX` taken from `fit_km`). |
|
1144 |
#' |
|
1145 |
#' @examples |
|
1146 |
#' \donttest{ |
|
1147 |
#' library(dplyr) |
|
1148 |
#' library(survival) |
|
1149 |
#' |
|
1150 |
#' adtte <- tern_ex_adtte %>% filter(PARAMCD == "OS") |
|
1151 |
#' fit <- survfit( |
|
1152 |
#' form = Surv(AVAL, 1 - CNSR) ~ ARMCD, |
|
1153 |
#' data = adtte |
|
1154 |
#' ) |
|
1155 |
#' h_tbl_median_surv(fit_km = fit) |
|
1156 |
#' } |
|
1157 |
#' |
|
1158 |
#' @export |
|
1159 |
h_tbl_median_surv <- function(fit_km, armval = "All") { |
|
1160 | 6x |
y <- if (is.null(fit_km$strata)) { |
1161 | ! |
as.data.frame(t(summary(fit_km)$table), row.names = armval) |
1162 |
} else { |
|
1163 | 6x |
tbl <- summary(fit_km)$table |
1164 | 6x |
rownames_lst <- strsplit(sub("=", "equals", rownames(tbl)), "equals") |
1165 | 6x |
rownames(tbl) <- matrix(unlist(rownames_lst), ncol = 2, byrow = TRUE)[, 2] |
1166 | 6x |
as.data.frame(tbl) |
1167 |
} |
|
1168 | 6x |
conf.int <- summary(fit_km)$conf.int # nolint |
1169 | 6x |
y$records <- round(y$records) |
1170 | 6x |
y$median <- signif(y$median, 4) |
1171 | 6x |
y$`CI` <- paste0( |
1172 | 6x |
"(", signif(y[[paste0(conf.int, "LCL")]], 4), ", ", signif(y[[paste0(conf.int, "UCL")]], 4), ")" |
1173 |
) |
|
1174 | 6x |
stats::setNames( |
1175 | 6x |
y[c("records", "median", "CI")], |
1176 | 6x |
c("N", "Median", f_conf_level(conf.int)) |
1177 |
) |
|
1178 |
} |
|
1179 | ||
1180 |
#' Helper Function: Survival Estimation Grob |
|
1181 |
#' |
|
1182 |
#' @description `r lifecycle::badge("stable")` |
|
1183 |
#' |
|
1184 |
#' The survival fit is transformed in a grob containing a table with groups in |
|
1185 |
#' rows characterized by N, median and 95% confidence interval. |
|
1186 |
#' |
|
1187 |
#' @inheritParams g_km |
|
1188 |
#' @inheritParams h_data_plot |
|
1189 |
#' @param ttheme (`list`)\cr see [gridExtra::ttheme_default()]. |
|
1190 |
#' @param x (`numeric`)\cr a value between 0 and 1 specifying x-location. |
|
1191 |
#' @param y (`numeric`)\cr a value between 0 and 1 specifying y-location. |
|
1192 |
#' @param width (`unit`)\cr width (as a unit) to use when printing the grob. |
|
1193 |
#' |
|
1194 |
#' @return A `grob` of a table containing statistics `N`, `Median`, and `XX% CI` (`XX` taken from `fit_km`). |
|
1195 |
#' |
|
1196 |
#' @examples |
|
1197 |
#' \donttest{ |
|
1198 |
#' library(dplyr) |
|
1199 |
#' library(survival) |
|
1200 |
#' library(grid) |
|
1201 |
#' |
|
1202 |
#' grid::grid.newpage() |
|
1203 |
#' grid.rect(gp = grid::gpar(lty = 1, col = "pink", fill = "gray85", lwd = 1)) |
|
1204 |
#' tern_ex_adtte %>% |
|
1205 |
#' filter(PARAMCD == "OS") %>% |
|
1206 |
#' survfit(form = Surv(AVAL, 1 - CNSR) ~ ARMCD, data = .) %>% |
|
1207 |
#' h_grob_median_surv() %>% |
|
1208 |
#' grid::grid.draw() |
|
1209 |
#' } |
|
1210 |
#' |
|
1211 |
#' @export |
|
1212 |
h_grob_median_surv <- function(fit_km, |
|
1213 |
armval = "All", |
|
1214 |
x = 0.9, |
|
1215 |
y = 0.9, |
|
1216 |
width = grid::unit(0.3, "npc"), |
|
1217 |
ttheme = gridExtra::ttheme_default()) { |
|
1218 | 5x |
data <- h_tbl_median_surv(fit_km, armval = armval) |
1219 | ||
1220 | 5x |
width <- grid::convertUnit(width, "in") |
1221 | 5x |
height <- width * (nrow(data) + 1) / 12 |
1222 | ||
1223 | 5x |
w <- paste(" ", c( |
1224 | 5x |
rownames(data)[which.max(nchar(rownames(data)))], |
1225 | 5x |
sapply(names(data), function(x) c(x, data[[x]])[which.max(nchar(c(x, data[[x]])))]) |
1226 |
)) |
|
1227 | 5x |
w_unit <- grid::convertWidth(grid::stringWidth(w), "in", valueOnly = TRUE) |
1228 | ||
1229 | 5x |
w_txt <- sapply(1:64, function(x) { |
1230 | 320x |
graphics::par(ps = x) |
1231 | 320x |
graphics::strwidth(w[4], units = "in") |
1232 |
}) |
|
1233 | 5x |
f_size_w <- which.max(w_txt[w_txt < as.numeric((w_unit / sum(w_unit)) * width)[4]]) |
1234 | ||
1235 | 5x |
h_txt <- sapply(1:64, function(x) { |
1236 | 320x |
graphics::par(ps = x) |
1237 | 320x |
graphics::strheight(grid::stringHeight("X"), units = "in") |
1238 |
}) |
|
1239 | 5x |
f_size_h <- which.max(h_txt[h_txt < as.numeric(grid::unit(as.numeric(height) / 4, grid::unitType(height)))]) |
1240 | ||
1241 | 5x |
if (ttheme$core$fg_params$fontsize == 12) { |
1242 | 5x |
ttheme$core$fg_params$fontsize <- min(f_size_w, f_size_h) |
1243 | 5x |
ttheme$colhead$fg_params$fontsize <- min(f_size_w, f_size_h) |
1244 | 5x |
ttheme$rowhead$fg_params$fontsize <- min(f_size_w, f_size_h) |
1245 |
} |
|
1246 | ||
1247 | 5x |
gt <- gridExtra::tableGrob( |
1248 | 5x |
d = data, |
1249 | 5x |
theme = ttheme |
1250 |
) |
|
1251 | 5x |
gt$widths <- ((w_unit / sum(w_unit)) * width) |
1252 | 5x |
gt$heights <- rep(grid::unit(as.numeric(height) / 4, grid::unitType(height)), nrow(gt)) |
1253 | ||
1254 | 5x |
vp <- grid::viewport( |
1255 | 5x |
x = grid::unit(x, "npc") + grid::unit(1, "lines"), |
1256 | 5x |
y = grid::unit(y, "npc") + grid::unit(1.5, "lines"), |
1257 | 5x |
height = height, |
1258 | 5x |
width = width, |
1259 | 5x |
just = c("right", "top") |
1260 |
) |
|
1261 | ||
1262 | 5x |
grid::gList( |
1263 | 5x |
grid::gTree( |
1264 | 5x |
vp = vp, |
1265 | 5x |
children = grid::gList(gt) |
1266 |
) |
|
1267 |
) |
|
1268 |
} |
|
1269 | ||
1270 |
#' Helper: Grid Object with y-axis Annotation |
|
1271 |
#' |
|
1272 |
#' @description `r lifecycle::badge("stable")` |
|
1273 |
#' |
|
1274 |
#' Build the y-axis annotation from a decomposed `ggplot`. |
|
1275 |
#' |
|
1276 |
#' @param ylab (`gtable`)\cr the y-lab as a graphical object derived from a `ggplot`. |
|
1277 |
#' @param yaxis (`gtable`)\cr the y-axis as a graphical object derived from a `ggplot`. |
|
1278 |
#' |
|
1279 |
#' @return a `gTree` object containing the y-axis annotation from a `ggplot`. |
|
1280 |
#' |
|
1281 |
#' @examples |
|
1282 |
#' \donttest{ |
|
1283 |
#' library(dplyr) |
|
1284 |
#' library(survival) |
|
1285 |
#' library(grid) |
|
1286 |
#' |
|
1287 |
#' fit_km <- tern_ex_adtte %>% |
|
1288 |
#' filter(PARAMCD == "OS") %>% |
|
1289 |
#' survfit(form = Surv(AVAL, 1 - CNSR) ~ ARMCD, data = .) |
|
1290 |
#' data_plot <- h_data_plot(fit_km = fit_km) |
|
1291 |
#' xticks <- h_xticks(data = data_plot) |
|
1292 |
#' gg <- h_ggkm( |
|
1293 |
#' data = data_plot, |
|
1294 |
#' censor_show = TRUE, |
|
1295 |
#' xticks = xticks, xlab = "Days", ylab = "Survival Probability", |
|
1296 |
#' title = "title", footnotes = "footnotes", yval = "Survival" |
|
1297 |
#' ) |
|
1298 |
#' |
|
1299 |
#' g_el <- h_decompose_gg(gg) |
|
1300 |
#' |
|
1301 |
#' grid::grid.newpage() |
|
1302 |
#' pvp <- grid::plotViewport(margins = c(5, 4, 2, 20)) |
|
1303 |
#' pushViewport(pvp) |
|
1304 |
#' grid::grid.draw(h_grob_y_annot(ylab = g_el$ylab, yaxis = g_el$yaxis)) |
|
1305 |
#' grid.rect(gp = grid::gpar(lty = 1, col = "gray35", fill = NA)) |
|
1306 |
#' } |
|
1307 |
#' |
|
1308 |
#' @export |
|
1309 |
h_grob_y_annot <- function(ylab, yaxis) { |
|
1310 | 6x |
grid::gList( |
1311 | 6x |
grid::gTree( |
1312 | 6x |
vp = grid::viewport( |
1313 | 6x |
width = grid::convertX(yaxis$width + ylab$width, "pt"), |
1314 | 6x |
x = grid::unit(1, "npc"), |
1315 | 6x |
just = "right" |
1316 |
), |
|
1317 | 6x |
children = grid::gList(cbind(ylab, yaxis)) |
1318 |
) |
|
1319 |
) |
|
1320 |
} |
|
1321 | ||
1322 |
#' Helper Function: Pairwise `CoxPH` table |
|
1323 |
#' |
|
1324 |
#' @description `r lifecycle::badge("stable")` |
|
1325 |
#' |
|
1326 |
#' Create a `data.frame` of pairwise stratified or unstratified `CoxPH` analysis results. |
|
1327 |
#' |
|
1328 |
#' @inheritParams g_km |
|
1329 |
#' |
|
1330 |
#' @return A `data.frame` containing statistics `HR`, `XX% CI` (`XX` taken from `control_coxph_pw`), |
|
1331 |
#' and `p-value (log-rank)`. |
|
1332 |
#' |
|
1333 |
#' @examples |
|
1334 |
#' \donttest{ |
|
1335 |
#' library(dplyr) |
|
1336 |
#' |
|
1337 |
#' adtte <- tern_ex_adtte %>% |
|
1338 |
#' filter(PARAMCD == "OS") %>% |
|
1339 |
#' mutate(is_event = CNSR == 0) |
|
1340 |
#' |
|
1341 |
#' h_tbl_coxph_pairwise( |
|
1342 |
#' df = adtte, |
|
1343 |
#' variables = list(tte = "AVAL", is_event = "is_event", arm = "ARM"), |
|
1344 |
#' control_coxph_pw = control_coxph(conf_level = 0.9) |
|
1345 |
#' ) |
|
1346 |
#' } |
|
1347 |
#' |
|
1348 |
#' @export |
|
1349 |
h_tbl_coxph_pairwise <- function(df, |
|
1350 |
variables, |
|
1351 |
control_coxph_pw = control_coxph()) { |
|
1352 | 3x |
assert_df_with_variables(df, variables) |
1353 | 3x |
arm <- variables$arm |
1354 | 3x |
df[[arm]] <- factor(df[[arm]]) |
1355 | 3x |
ref_group <- levels(df[[arm]])[1] |
1356 | 3x |
comp_group <- levels(df[[arm]])[-1] |
1357 | 3x |
results <- Map(function(comp) { |
1358 | 6x |
res <- s_coxph_pairwise( |
1359 | 6x |
df = df[df[[arm]] == comp, , drop = FALSE], |
1360 | 6x |
.ref_group = df[df[[arm]] == ref_group, , drop = FALSE], |
1361 | 6x |
.in_ref_col = FALSE, |
1362 | 6x |
.var = variables$tte, |
1363 | 6x |
is_event = variables$is_event, |
1364 | 6x |
strat = variables$strat, |
1365 | 6x |
control = control_coxph_pw |
1366 |
) |
|
1367 | 6x |
res_df <- data.frame( |
1368 | 6x |
hr = format(round(res$hr, 2), nsmall = 2), |
1369 | 6x |
hr_ci = paste0( |
1370 | 6x |
"(", format(round(res$hr_ci[1], 2), nsmall = 2), ", ", |
1371 | 6x |
format(round(res$hr_ci[2], 2), nsmall = 2), ")" |
1372 |
), |
|
1373 | 6x |
pvalue = if (res$pvalue < 0.0001) "<0.0001" else format(round(res$pvalue, 4), 4), |
1374 | 6x |
stringsAsFactors = FALSE |
1375 |
) |
|
1376 | 6x |
colnames(res_df) <- c("HR", vapply(res[c("hr_ci", "pvalue")], obj_label, FUN.VALUE = "character")) |
1377 | 6x |
row.names(res_df) <- comp |
1378 | 6x |
res_df |
1379 | 3x |
}, comp_group) |
1380 | 3x |
do.call(rbind, results) |
1381 |
} |
|
1382 | ||
1383 |
#' Helper Function: `CoxPH` Grob |
|
1384 |
#' |
|
1385 |
#' @description `r lifecycle::badge("stable")` |
|
1386 |
#' |
|
1387 |
#' Grob of `rtable` output from [h_tbl_coxph_pairwise()] |
|
1388 |
#' |
|
1389 |
#' @inheritParams h_grob_median_surv |
|
1390 |
#' @param ... arguments will be passed to [h_tbl_coxph_pairwise()]. |
|
1391 |
#' @param x (`numeric`)\cr a value between 0 and 1 specifying x-location. |
|
1392 |
#' @param y (`numeric`)\cr a value between 0 and 1 specifying y-location. |
|
1393 |
#' @param width (`unit`)\cr width (as a unit) to use when printing the grob. |
|
1394 |
#' |
|
1395 |
#' @return A `grob` of a table containing statistics `HR`, `XX% CI` (`XX` taken from `control_coxph_pw`), |
|
1396 |
#' and `p-value (log-rank)`. |
|
1397 |
#' |
|
1398 |
#' @examples |
|
1399 |
#' \donttest{ |
|
1400 |
#' library(dplyr) |
|
1401 |
#' library(survival) |
|
1402 |
#' library(grid) |
|
1403 |
#' |
|
1404 |
#' grid::grid.newpage() |
|
1405 |
#' grid.rect(gp = grid::gpar(lty = 1, col = "pink", fill = "gray85", lwd = 1)) |
|
1406 |
#' data <- tern_ex_adtte %>% |
|
1407 |
#' filter(PARAMCD == "OS") %>% |
|
1408 |
#' mutate(is_event = CNSR == 0) |
|
1409 |
#' tbl_grob <- h_grob_coxph( |
|
1410 |
#' df = data, |
|
1411 |
#' variables = list(tte = "AVAL", is_event = "is_event", arm = "ARMCD"), |
|
1412 |
#' control_coxph_pw = control_coxph(conf_level = 0.9), x = 0.5, y = 0.5 |
|
1413 |
#' ) |
|
1414 |
#' grid::grid.draw(tbl_grob) |
|
1415 |
#' } |
|
1416 |
#' |
|
1417 |
#' @export |
|
1418 |
h_grob_coxph <- function(..., |
|
1419 |
x = 0, |
|
1420 |
y = 0, |
|
1421 |
width = grid::unit(0.4, "npc"), |
|
1422 |
ttheme = gridExtra::ttheme_default( |
|
1423 |
padding = grid::unit(c(1, .5), "lines"), |
|
1424 |
core = list(bg_params = list(fill = c("grey95", "grey90"), alpha = .5)) |
|
1425 |
)) { |
|
1426 | 2x |
data <- h_tbl_coxph_pairwise(...) |
1427 | ||
1428 | 2x |
width <- grid::convertUnit(width, "in") |
1429 | 2x |
height <- width * (nrow(data) + 1) / 12 |
1430 | ||
1431 | 2x |
w <- paste(" ", c( |
1432 | 2x |
rownames(data)[which.max(nchar(rownames(data)))], |
1433 | 2x |
sapply(names(data), function(x) c(x, data[[x]])[which.max(nchar(c(x, data[[x]])))]) |
1434 |
)) |
|
1435 | 2x |
w_unit <- grid::convertWidth(grid::stringWidth(w), "in", valueOnly = TRUE) |
1436 | ||
1437 | 2x |
w_txt <- sapply(1:64, function(x) { |
1438 | 128x |
graphics::par(ps = x) |
1439 | 128x |
graphics::strwidth(w[4], units = "in") |
1440 |
}) |
|
1441 | 2x |
f_size_w <- which.max(w_txt[w_txt < as.numeric((w_unit / sum(w_unit)) * width)[4]]) |
1442 | ||
1443 | 2x |
h_txt <- sapply(1:64, function(x) { |
1444 | 128x |
graphics::par(ps = x) |
1445 | 128x |
graphics::strheight(grid::stringHeight("X"), units = "in") |
1446 |
}) |
|
1447 | 2x |
f_size_h <- which.max(h_txt[h_txt < as.numeric(grid::unit(as.numeric(height) / 4, grid::unitType(height)))]) |
1448 | ||
1449 | 2x |
if (ttheme$core$fg_params$fontsize == 12) { |
1450 | 2x |
ttheme$core$fg_params$fontsize <- min(f_size_w, f_size_h) |
1451 | 2x |
ttheme$colhead$fg_params$fontsize <- min(f_size_w, f_size_h) |
1452 | 2x |
ttheme$rowhead$fg_params$fontsize <- min(f_size_w, f_size_h) |
1453 |
} |
|
1454 | ||
1455 | 2x |
tryCatch( |
1456 | 2x |
expr = { |
1457 | 2x |
gt <- gridExtra::tableGrob( |
1458 | 2x |
d = data, |
1459 | 2x |
theme = ttheme |
1460 | 2x |
) # ERROR 'data' must be of a vector type, was 'NULL' |
1461 | 2x |
gt$widths <- ((w_unit / sum(w_unit)) * width) |
1462 | 2x |
gt$heights <- rep(grid::unit(as.numeric(height) / 4, grid::unitType(height)), nrow(gt)) |
1463 | 2x |
vp <- grid::viewport( |
1464 | 2x |
x = grid::unit(x, "npc") + grid::unit(1, "lines"), |
1465 | 2x |
y = grid::unit(y, "npc") + grid::unit(1.5, "lines"), |
1466 | 2x |
height = height, |
1467 | 2x |
width = width, |
1468 | 2x |
just = c("left", "bottom") |
1469 |
) |
|
1470 | 2x |
grid::gList( |
1471 | 2x |
grid::gTree( |
1472 | 2x |
vp = vp, |
1473 | 2x |
children = grid::gList(gt) |
1474 |
) |
|
1475 |
) |
|
1476 |
}, |
|
1477 | 2x |
error = function(w) { |
1478 | ! |
message(paste( |
1479 | ! |
"Warning: Cox table will not be displayed as there is", |
1480 | ! |
"not any level to be compared in the arm variable." |
1481 |
)) |
|
1482 | ! |
return( |
1483 | ! |
grid::gList( |
1484 | ! |
grid::gTree( |
1485 | ! |
vp = NULL, |
1486 | ! |
children = NULL |
1487 |
) |
|
1488 |
) |
|
1489 |
) |
|
1490 |
} |
|
1491 |
) |
|
1492 |
} |
1 |
#' Add Titles, Footnotes, Page Number, and a Bounding Box to a Grid Grob |
|
2 |
#' |
|
3 |
#' @description `r lifecycle::badge("stable")` |
|
4 |
#' |
|
5 |
#' This function is useful to label grid grobs (also `ggplot2`, and `lattice` plots) |
|
6 |
#' with title, footnote, and page numbers. |
|
7 |
#' |
|
8 |
#' @inheritParams grid::grob |
|
9 |
#' @param grob a grid grob object, optionally `NULL` if only a `grob` with the decoration should be shown. |
|
10 |
#' @param titles vector of character strings. Vector elements are separated by a newline and strings are wrapped |
|
11 |
#' according to the page width. |
|
12 |
#' @param footnotes vector of character string. Same rules as for `titles`. |
|
13 |
#' @param page string with page numeration, if `NULL` then no page number is displayed. |
|
14 |
#' @param width_titles unit object |
|
15 |
#' @param width_footnotes unit object |
|
16 |
#' @param border boolean, whether a a border should be drawn around the plot or not. |
|
17 |
#' @param margins unit object of length 4 |
|
18 |
#' @param padding unit object of length 4 |
|
19 |
#' @param outer_margins unit object of length 4 |
|
20 |
#' @param gp_titles a `gpar` object |
|
21 |
#' @param gp_footnotes a `gpar` object |
|
22 |
#' |
|
23 |
#' @return A grid grob (`gTree`). |
|
24 |
#' |
|
25 |
#' @details The titles and footnotes will be ragged, i.e. each title will be wrapped individually. |
|
26 |
#' |
|
27 |
#' @examples |
|
28 |
#' library(grid) |
|
29 |
#' |
|
30 |
#' titles <- c( |
|
31 |
#' "Edgar Anderson's Iris Data", |
|
32 |
#' paste( |
|
33 |
#' "This famous (Fisher's or Anderson's) iris data set gives the measurements", |
|
34 |
#' "in centimeters of the variables sepal length and width and petal length", |
|
35 |
#' "and width, respectively, for 50 flowers from each of 3 species of iris." |
|
36 |
#' ) |
|
37 |
#' ) |
|
38 |
#' |
|
39 |
#' footnotes <- c( |
|
40 |
#' "The species are Iris setosa, versicolor, and virginica.", |
|
41 |
#' paste( |
|
42 |
#' "iris is a data frame with 150 cases (rows) and 5 variables (columns) named", |
|
43 |
#' "Sepal.Length, Sepal.Width, Petal.Length, Petal.Width, and Species." |
|
44 |
#' ) |
|
45 |
#' ) |
|
46 |
#' |
|
47 |
#' ## empty plot |
|
48 |
#' grid.newpage() |
|
49 |
#' |
|
50 |
#' grid.draw( |
|
51 |
#' decorate_grob( |
|
52 |
#' NULL, |
|
53 |
#' titles = titles, |
|
54 |
#' footnotes = footnotes, |
|
55 |
#' page = "Page 4 of 10" |
|
56 |
#' ) |
|
57 |
#' ) |
|
58 |
#' |
|
59 |
#' # grid |
|
60 |
#' p <- gTree( |
|
61 |
#' children = gList( |
|
62 |
#' rectGrob(), |
|
63 |
#' xaxisGrob(), |
|
64 |
#' yaxisGrob(), |
|
65 |
#' textGrob("Sepal.Length", y = unit(-4, "lines")), |
|
66 |
#' textGrob("Petal.Length", x = unit(-3.5, "lines"), rot = 90), |
|
67 |
#' pointsGrob(iris$Sepal.Length, iris$Petal.Length, gp = gpar(col = iris$Species), pch = 16) |
|
68 |
#' ), |
|
69 |
#' vp = vpStack(plotViewport(), dataViewport(xData = iris$Sepal.Length, yData = iris$Petal.Length)) |
|
70 |
#' ) |
|
71 |
#' grid.newpage() |
|
72 |
#' grid.draw(p) |
|
73 |
#' |
|
74 |
#' grid.newpage() |
|
75 |
#' grid.draw( |
|
76 |
#' decorate_grob( |
|
77 |
#' grob = p, |
|
78 |
#' titles = titles, |
|
79 |
#' footnotes = footnotes, |
|
80 |
#' page = "Page 6 of 129" |
|
81 |
#' ) |
|
82 |
#' ) |
|
83 |
#' |
|
84 |
#' ## with ggplot2 |
|
85 |
#' library(ggplot2) |
|
86 |
#' |
|
87 |
#' p_gg <- ggplot2::ggplot(iris, aes(Sepal.Length, Sepal.Width, col = Species)) + |
|
88 |
#' ggplot2::geom_point() |
|
89 |
#' p_gg |
|
90 |
#' p <- ggplotGrob(p_gg) |
|
91 |
#' grid.newpage() |
|
92 |
#' grid.draw( |
|
93 |
#' decorate_grob( |
|
94 |
#' grob = p, |
|
95 |
#' titles = titles, |
|
96 |
#' footnotes = footnotes, |
|
97 |
#' page = "Page 6 of 129" |
|
98 |
#' ) |
|
99 |
#' ) |
|
100 |
#' |
|
101 |
#' ## with lattice |
|
102 |
#' library(lattice) |
|
103 |
#' |
|
104 |
#' xyplot(Sepal.Length ~ Petal.Length, data = iris, col = iris$Species) |
|
105 |
#' p <- grid.grab() |
|
106 |
#' grid.newpage() |
|
107 |
#' grid.draw( |
|
108 |
#' decorate_grob( |
|
109 |
#' grob = p, |
|
110 |
#' titles = titles, |
|
111 |
#' footnotes = footnotes, |
|
112 |
#' page = "Page 6 of 129" |
|
113 |
#' ) |
|
114 |
#' ) |
|
115 |
#' |
|
116 |
#' # with gridExtra - no borders |
|
117 |
#' library(gridExtra) |
|
118 |
#' grid.newpage() |
|
119 |
#' grid.draw( |
|
120 |
#' decorate_grob( |
|
121 |
#' tableGrob( |
|
122 |
#' head(mtcars) |
|
123 |
#' ), |
|
124 |
#' titles = "title", |
|
125 |
#' footnotes = "footnote", |
|
126 |
#' border = FALSE |
|
127 |
#' ) |
|
128 |
#' ) |
|
129 |
#' |
|
130 |
#' @export |
|
131 |
decorate_grob <- function(grob, |
|
132 |
titles, |
|
133 |
footnotes, |
|
134 |
page = "", |
|
135 |
width_titles = grid::unit(1, "npc") - grid::stringWidth(page), |
|
136 |
width_footnotes = grid::unit(1, "npc") - grid::stringWidth(page), |
|
137 |
border = TRUE, |
|
138 |
margins = grid::unit(c(1, 0, 1, 0), "lines"), |
|
139 |
padding = grid::unit(rep(1, 4), "lines"), |
|
140 |
outer_margins = grid::unit(c(2, 1.5, 3, 1.5), "cm"), |
|
141 |
gp_titles = grid::gpar(), |
|
142 |
gp_footnotes = grid::gpar(fontsize = 8), |
|
143 |
name = NULL, |
|
144 |
gp = grid::gpar(), |
|
145 |
vp = NULL) { |
|
146 | 8x |
st_titles <- split_text_grob( |
147 | 8x |
titles, |
148 | 8x |
x = 0, y = 1, |
149 | 8x |
just = c("left", "top"), |
150 | 8x |
width = width_titles, |
151 | 8x |
vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 1), |
152 | 8x |
gp = gp_titles |
153 |
) |
|
154 | ||
155 | 8x |
st_footnotes <- split_text_grob( |
156 | 8x |
footnotes, |
157 | 8x |
x = 0, y = 1, |
158 | 8x |
just = c("left", "top"), |
159 | 8x |
width = width_footnotes, |
160 | 8x |
vp = grid::viewport(layout.pos.row = 3, layout.pos.col = 1), |
161 | 8x |
gp = gp_footnotes |
162 |
) |
|
163 | ||
164 | 8x |
grid::gTree( |
165 | 8x |
grob = grob, |
166 | 8x |
titles = titles, |
167 | 8x |
footnotes = footnotes, |
168 | 8x |
page = page, |
169 | 8x |
width_titles = width_titles, |
170 | 8x |
width_footnotes = width_footnotes, |
171 | 8x |
border = border, |
172 | 8x |
margins = margins, |
173 | 8x |
padding = padding, |
174 | 8x |
outer_margins = outer_margins, |
175 | 8x |
gp_titles = gp_titles, |
176 | 8x |
gp_footnotes = gp_footnotes, |
177 | 8x |
children = grid::gList( |
178 | 8x |
grid::gTree( |
179 | 8x |
children = grid::gList( |
180 | 8x |
st_titles, |
181 | 8x |
grid::gTree( |
182 | 8x |
children = grid::gList( |
183 | 8x |
if (border) grid::rectGrob(), |
184 | 8x |
grid::gTree( |
185 | 8x |
children = grid::gList( |
186 | 8x |
grob |
187 |
), |
|
188 | 8x |
vp = grid::plotViewport(margins = padding) |
189 |
) |
|
190 |
), |
|
191 | 8x |
vp = grid::vpStack( |
192 | 8x |
grid::viewport(layout.pos.row = 2, layout.pos.col = 1), |
193 | 8x |
grid::plotViewport(margins = margins) |
194 |
) |
|
195 |
), |
|
196 | 8x |
st_footnotes, |
197 | 8x |
grid::textGrob( |
198 | 8x |
page, |
199 | 8x |
x = 1, y = 0, |
200 | 8x |
just = c("right", "bottom"), |
201 | 8x |
vp = grid::viewport(layout.pos.row = 3, layout.pos.col = 1), |
202 | 8x |
gp = gp_footnotes |
203 |
) |
|
204 |
), |
|
205 | 8x |
childrenvp = NULL, |
206 | 8x |
name = "titles_grob_footnotes", |
207 | 8x |
vp = grid::vpStack( |
208 | 8x |
grid::plotViewport(margins = outer_margins), |
209 | 8x |
grid::viewport( |
210 | 8x |
layout = grid::grid.layout( |
211 | 8x |
nrow = 3, ncol = 1, |
212 | 8x |
heights = grid::unit.c( |
213 | 8x |
grid::grobHeight(st_titles), |
214 | 8x |
grid::unit(1, "null"), |
215 | 8x |
grid::grobHeight(st_footnotes) |
216 |
) |
|
217 |
) |
|
218 |
) |
|
219 |
) |
|
220 |
) |
|
221 |
), |
|
222 | 8x |
name = name, |
223 | 8x |
gp = gp, |
224 | 8x |
vp = vp, |
225 | 8x |
cl = "decoratedGrob" |
226 |
) |
|
227 |
} |
|
228 | ||
229 |
#' @importFrom grid validDetails |
|
230 |
#' @noRd |
|
231 |
validDetails.decoratedGrob <- function(x) { |
|
232 | ! |
checkmate::assert_character(x$titles) |
233 | ! |
checkmate::assert_character(x$footnotes) |
234 | ||
235 | ! |
if (!is.null(x$grob)) { |
236 | ! |
checkmate::assert_true(grid::is.grob(x$grob)) |
237 |
} |
|
238 | ! |
if (length(x$page) == 1) { |
239 | ! |
checkmate::assert_character(x$page) |
240 |
} |
|
241 | ! |
if (!grid::is.unit(x$outer_margins)) { |
242 | ! |
checkmate::assert_vector(x$outer_margins, len = 4) |
243 |
} |
|
244 | ! |
if (!grid::is.unit(x$margins)) { |
245 | ! |
checkmate::assert_vector(x$margins, len = 4) |
246 |
} |
|
247 | ! |
if (!grid::is.unit(x$padding)) { |
248 | ! |
checkmate::assert_vector(x$padding, len = 4) |
249 |
} |
|
250 | ||
251 | ! |
x |
252 |
} |
|
253 | ||
254 |
#' @importFrom grid widthDetails |
|
255 |
#' @noRd |
|
256 |
widthDetails.decoratedGrob <- function(x) { |
|
257 | ! |
grid::unit(1, "null") |
258 |
} |
|
259 | ||
260 |
#' @importFrom grid heightDetails |
|
261 |
#' @noRd |
|
262 |
heightDetails.decoratedGrob <- function(x) { |
|
263 | ! |
grid::unit(1, "null") |
264 |
} |
|
265 | ||
266 |
# Adapted from Paul Murell R Graphics 2nd Edition |
|
267 |
# https://www.stat.auckland.ac.nz/~paul/RG2e/interactgrid-splittext.R |
|
268 |
split_string <- function(text, width) { |
|
269 | 17x |
strings <- strsplit(text, " ") |
270 | 17x |
out_string <- NA |
271 | 17x |
for (string_i in seq_along(strings)) { |
272 | 17x |
newline_str <- strings[[string_i]] |
273 | 6x |
if (length(newline_str) == 0) newline_str <- "" |
274 | 17x |
if (is.na(out_string[string_i])) { |
275 | 17x |
out_string[string_i] <- newline_str[[1]][[1]] |
276 | 17x |
linewidth <- grid::stringWidth(out_string[string_i]) |
277 |
} |
|
278 | 17x |
gapwidth <- grid::stringWidth(" ") |
279 | 17x |
availwidth <- as.numeric(width) |
280 | 17x |
if (length(newline_str) > 1) { |
281 | 5x |
for (i in seq(2, length(newline_str))) { |
282 | 27x |
width_i <- grid::stringWidth(newline_str[i]) |
283 | 27x |
if (grid::convertWidth(linewidth + gapwidth + width_i, grid::unitType(width), valueOnly = TRUE) < availwidth) { |
284 | 25x |
sep <- " " |
285 | 25x |
linewidth <- linewidth + gapwidth + width_i |
286 |
} else { |
|
287 | 2x |
sep <- "\n" |
288 | 2x |
linewidth <- width_i |
289 |
} |
|
290 | 27x |
out_string[string_i] <- paste(out_string[string_i], newline_str[i], sep = sep) |
291 |
} |
|
292 |
} |
|
293 |
} |
|
294 | 17x |
paste(out_string, collapse = "\n") |
295 |
} |
|
296 | ||
297 |
#' Split Text According To Available Text Width |
|
298 |
#' |
|
299 |
#' Dynamically wrap text. |
|
300 |
#' |
|
301 |
#' @inheritParams grid::grid.text |
|
302 |
#' @param text character string |
|
303 |
#' @param width a unit object specifying max width of text |
|
304 |
#' |
|
305 |
#' @return A text grob. |
|
306 |
#' |
|
307 |
#' @details This code is taken from `R Graphics by Paul Murell, 2nd edition` |
|
308 |
#' |
|
309 |
#' @keywords internal |
|
310 |
split_text_grob <- function(text, |
|
311 |
x = grid::unit(0.5, "npc"), |
|
312 |
y = grid::unit(0.5, "npc"), |
|
313 |
width = grid::unit(1, "npc"), |
|
314 |
just = "centre", |
|
315 |
hjust = NULL, |
|
316 |
vjust = NULL, |
|
317 |
default.units = "npc", # nolint |
|
318 |
name = NULL, |
|
319 |
gp = grid::gpar(), |
|
320 |
vp = NULL) { |
|
321 | 16x |
if (!grid::is.unit(x)) x <- grid::unit(x, default.units) |
322 | 16x |
if (!grid::is.unit(y)) y <- grid::unit(y, default.units) |
323 | ! |
if (!grid::is.unit(width)) width <- grid::unit(width, default.units) |
324 | ! |
if (grid::unitType(x) %in% c("sum", "min", "max")) x <- grid::convertUnit(x, default.units) |
325 | ! |
if (grid::unitType(y) %in% c("sum", "min", "max")) y <- grid::convertUnit(y, default.units) |
326 | 16x |
if (grid::unitType(width) %in% c("sum", "min", "max")) width <- grid::convertUnit(width, default.units) |
327 | ||
328 |
## if it is a fixed unit then we do not need to recalculate when viewport resized |
|
329 | 16x |
if (!inherits(width, "unit.arithmetic") && |
330 | 16x |
!is.null(attr(width, "unit")) && |
331 | 16x |
attr(width, "unit") %in% c("cm", "inches", "mm", "points", "picas", "bigpts", "dida", "cicero", "scaledpts")) { |
332 | ! |
attr(text, "fixed_text") <- paste(vapply(text, split_string, character(1), width = width), collapse = "\n") |
333 |
} |
|
334 | ||
335 | 16x |
grid::grid.text( |
336 | 16x |
label = split_string(text, width), |
337 | 16x |
x = x, y = y, |
338 | 16x |
just = just, |
339 | 16x |
hjust = hjust, |
340 | 16x |
vjust = vjust, |
341 | 16x |
rot = 0, |
342 | 16x |
check.overlap = FALSE, |
343 | 16x |
name = name, |
344 | 16x |
gp = gp, |
345 | 16x |
vp = vp, |
346 | 16x |
draw = FALSE |
347 |
) |
|
348 |
} |
|
349 | ||
350 |
#' @importFrom grid validDetails |
|
351 |
#' @noRd |
|
352 |
validDetails.dynamicSplitText <- function(x) { |
|
353 | ! |
checkmate::assert_character(x$text) |
354 | ! |
checkmate::assert_true(grid::is.unit(x$width)) |
355 | ! |
checkmate::assert_vector(x$width, len = 1) |
356 | ! |
x |
357 |
} |
|
358 | ||
359 |
#' @importFrom grid heightDetails |
|
360 |
#' @noRd |
|
361 |
heightDetails.dynamicSplitText <- function(x) { |
|
362 | ! |
txt <- if (!is.null(attr(x$text, "fixed_text"))) { |
363 | ! |
attr(x$text, "fixed_text") |
364 |
} else { |
|
365 | ! |
paste(vapply(x$text, split_string, character(1), width = x$width), collapse = "\n") |
366 |
} |
|
367 | ! |
grid::stringHeight(txt) |
368 |
} |
|
369 | ||
370 |
#' @importFrom grid widthDetails |
|
371 |
#' @noRd |
|
372 |
widthDetails.dynamicSplitText <- function(x) { |
|
373 | ! |
x$width |
374 |
} |
|
375 | ||
376 |
#' @importFrom grid drawDetails |
|
377 |
#' @noRd |
|
378 |
drawDetails.dynamicSplitText <- function(x, recording) { |
|
379 | ! |
txt <- if (!is.null(attr(x$text, "fixed_text"))) { |
380 | ! |
attr(x$text, "fixed_text") |
381 |
} else { |
|
382 | ! |
paste(vapply(x$text, split_string, character(1), width = x$width), collapse = "\n") |
383 |
} |
|
384 | ||
385 | ! |
x$width <- NULL |
386 | ! |
x$label <- txt |
387 | ! |
x$text <- NULL |
388 | ! |
class(x) <- c("text", class(x)[-1]) |
389 | ||
390 | ! |
grid::grid.draw(x) |
391 |
} |
|
392 | ||
393 |
#' Update Page Number |
|
394 |
#' |
|
395 |
#' Automatically updates page number. |
|
396 |
#' |
|
397 |
#' @param npages number of pages in total |
|
398 |
#' @param ... passed on to [decorate_grob()] |
|
399 |
#' |
|
400 |
#' @return Closure that increments the page number. |
|
401 |
#' |
|
402 |
#' @keywords internal |
|
403 |
decorate_grob_factory <- function(npages, ...) { |
|
404 | 2x |
current_page <- 0 |
405 | 2x |
function(grob) { |
406 | 7x |
current_page <<- current_page + 1 |
407 | 7x |
if (current_page > npages) { |
408 | 1x |
stop(paste("current page is", current_page, "but max.", npages, "specified.")) |
409 |
} |
|
410 | 6x |
decorate_grob(grob = grob, page = paste("Page", current_page, "of", npages), ...) |
411 |
} |
|
412 |
} |
|
413 | ||
414 |
#' Decorate Set of `grobs` and Add Page Numbering |
|
415 |
#' |
|
416 |
#' @description `r lifecycle::badge("stable")` |
|
417 |
#' |
|
418 |
#' Note that this uses the [decorate_grob_factory()] function. |
|
419 |
#' |
|
420 |
#' @param grobs a list of grid grobs |
|
421 |
#' @param ... arguments passed on to [decorate_grob()]. |
|
422 |
#' |
|
423 |
#' @return A decorated grob. |
|
424 |
#' |
|
425 |
#' @examples |
|
426 |
#' library(ggplot2) |
|
427 |
#' library(grid) |
|
428 |
#' g <- with(data = iris, { |
|
429 |
#' list( |
|
430 |
#' ggplot2::ggplotGrob( |
|
431 |
#' ggplot2::ggplot(mapping = aes(Sepal.Length, Sepal.Width, col = Species)) + |
|
432 |
#' ggplot2::geom_point() |
|
433 |
#' ), |
|
434 |
#' ggplot2::ggplotGrob( |
|
435 |
#' ggplot2::ggplot(mapping = aes(Sepal.Length, Petal.Length, col = Species)) + |
|
436 |
#' ggplot2::geom_point() |
|
437 |
#' ), |
|
438 |
#' ggplot2::ggplotGrob( |
|
439 |
#' ggplot2::ggplot(mapping = aes(Sepal.Length, Petal.Width, col = Species)) + |
|
440 |
#' ggplot2::geom_point() |
|
441 |
#' ), |
|
442 |
#' ggplot2::ggplotGrob( |
|
443 |
#' ggplot2::ggplot(mapping = aes(Sepal.Width, Petal.Length, col = Species)) + |
|
444 |
#' ggplot2::geom_point() |
|
445 |
#' ), |
|
446 |
#' ggplot2::ggplotGrob( |
|
447 |
#' ggplot2::ggplot(mapping = aes(Sepal.Width, Petal.Width, col = Species)) + |
|
448 |
#' ggplot2::geom_point() |
|
449 |
#' ), |
|
450 |
#' ggplot2::ggplotGrob( |
|
451 |
#' ggplot2::ggplot(mapping = aes(Petal.Length, Petal.Width, col = Species)) + |
|
452 |
#' ggplot2::geom_point() |
|
453 |
#' ) |
|
454 |
#' ) |
|
455 |
#' }) |
|
456 |
#' lg <- decorate_grob_set(grobs = g, titles = "Hello\nOne\nTwo\nThree", footnotes = "") |
|
457 |
#' |
|
458 |
#' draw_grob(lg[[1]]) |
|
459 |
#' draw_grob(lg[[2]]) |
|
460 |
#' draw_grob(lg[[6]]) |
|
461 |
#' |
|
462 |
#' @export |
|
463 |
decorate_grob_set <- function(grobs, ...) { |
|
464 | 1x |
n <- length(grobs) |
465 | 1x |
lgf <- decorate_grob_factory(npages = n, ...) |
466 | 1x |
lapply(grobs, lgf) |
467 |
} |
1 |
#' Confidence Intervals for a Difference of Binomials |
|
2 |
#' |
|
3 |
#' @description `r lifecycle::badge("experimental")` |
|
4 |
#' |
|
5 |
#' Several confidence intervals for the difference between proportions. |
|
6 |
#' |
|
7 |
#' @name desctools_binom |
|
8 |
NULL |
|
9 | ||
10 |
#' Recycle List of Parameters |
|
11 |
#' |
|
12 |
#' This function recycles all supplied elements to the maximal dimension. |
|
13 |
#' |
|
14 |
#' @param ... (`any`)\cr Elements to recycle. |
|
15 |
#' |
|
16 |
#' @return A `list`. |
|
17 |
#' |
|
18 |
#' @keywords internal |
|
19 |
#' @noRd |
|
20 |
h_recycle <- function(...) { |
|
21 | 60x |
lst <- list(...) |
22 | 60x |
maxdim <- max(lengths(lst)) |
23 | 60x |
res <- lapply(lst, rep, length.out = maxdim) |
24 | 60x |
attr(res, "maxdim") <- maxdim |
25 | 60x |
return(res) |
26 |
} |
|
27 | ||
28 |
#' @describeIn desctools_binom Several confidence intervals for the difference between proportions. |
|
29 |
#' |
|
30 |
#' @return A `matrix` of 3 values: |
|
31 |
#' * `est`: estimate of proportion difference. |
|
32 |
#' * `lwr.ci`: estimate of lower end of the confidence interval. |
|
33 |
#' * `upr.ci`: estimate of upper end of the confidence interval. |
|
34 |
#' |
|
35 |
#' @keywords internal |
|
36 |
desctools_binom <- function(x1, n1, x2, n2, conf.level = 0.95, sides = c( # nolint |
|
37 |
"two.sided", |
|
38 |
"left", "right" |
|
39 |
), method = c( |
|
40 |
"ac", "wald", "waldcc", "score", |
|
41 |
"scorecc", "mn", "mee", "blj", "ha", "hal", "jp" |
|
42 |
)) { |
|
43 | 18x |
if (missing(sides)) { |
44 | 18x |
sides <- match.arg(sides) |
45 |
} |
|
46 | 18x |
if (missing(method)) { |
47 | 1x |
method <- match.arg(method) |
48 |
} |
|
49 | 18x |
iBinomDiffCI <- function(x1, n1, x2, n2, conf.level, sides, # nolint |
50 | 18x |
method) { |
51 | 18x |
if (sides != "two.sided") { |
52 | ! |
conf.level <- 1 - 2 * (1 - conf.level) # nolint |
53 |
} |
|
54 | 18x |
alpha <- 1 - conf.level |
55 | 18x |
kappa <- stats::qnorm(1 - alpha / 2) |
56 | 18x |
p1_hat <- x1 / n1 |
57 | 18x |
p2_hat <- x2 / n2 |
58 | 18x |
est <- p1_hat - p2_hat |
59 | 18x |
switch(method, |
60 | 18x |
wald = { |
61 | 2x |
vd <- p1_hat * (1 - p1_hat) / n1 + p2_hat * (1 - p2_hat) / n2 |
62 | 2x |
term2 <- kappa * sqrt(vd) |
63 | 2x |
ci_lwr <- max(-1, est - term2) |
64 | 2x |
ci_upr <- min(1, est + term2) |
65 |
}, |
|
66 | 18x |
waldcc = { |
67 | 2x |
vd <- p1_hat * (1 - p1_hat) / n1 + p2_hat * (1 - p2_hat) / n2 |
68 | 2x |
term2 <- kappa * sqrt(vd) |
69 | 2x |
term2 <- term2 + 0.5 * (1 / n1 + 1 / n2) |
70 | 2x |
ci_lwr <- max(-1, est - term2) |
71 | 2x |
ci_upr <- min(1, est + term2) |
72 |
}, |
|
73 | 18x |
ac = { |
74 | 2x |
n1 <- n1 + 2 |
75 | 2x |
n2 <- n2 + 2 |
76 | 2x |
x1 <- x1 + 1 |
77 | 2x |
x2 <- x2 + 1 |
78 | 2x |
p1_hat <- x1 / n1 |
79 | 2x |
p2_hat <- x2 / n2 |
80 | 2x |
est1 <- p1_hat - p2_hat |
81 | 2x |
vd <- p1_hat * (1 - p1_hat) / n1 + p2_hat * (1 - p2_hat) / n2 |
82 | 2x |
term2 <- kappa * sqrt(vd) |
83 | 2x |
ci_lwr <- max(-1, est1 - term2) |
84 | 2x |
ci_upr <- min(1, est1 + term2) |
85 |
}, |
|
86 | 18x |
exact = { |
87 | ! |
ci_lwr <- NA |
88 | ! |
ci_upr <- NA |
89 |
}, |
|
90 | 18x |
score = { |
91 | 2x |
w1 <- desctools_binomci( |
92 | 2x |
x = x1, n = n1, conf.level = conf.level, |
93 | 2x |
method = "wilson" |
94 |
) |
|
95 | 2x |
w2 <- desctools_binomci( |
96 | 2x |
x = x2, n = n2, conf.level = conf.level, |
97 | 2x |
method = "wilson" |
98 |
) |
|
99 | 2x |
l1 <- w1[2] |
100 | 2x |
u1 <- w1[3] |
101 | 2x |
l2 <- w2[2] |
102 | 2x |
u2 <- w2[3] |
103 | 2x |
ci_lwr <- est - kappa * sqrt(l1 * (1 - l1) / n1 + |
104 | 2x |
u2 * (1 - u2) / n2) |
105 | 2x |
ci_upr <- est + kappa * sqrt(u1 * (1 - u1) / n1 + |
106 | 2x |
l2 * (1 - l2) / n2) |
107 |
}, |
|
108 | 18x |
scorecc = { |
109 | 1x |
w1 <- desctools_binomci( |
110 | 1x |
x = x1, n = n1, conf.level = conf.level, |
111 | 1x |
method = "wilsoncc" |
112 |
) |
|
113 | 1x |
w2 <- desctools_binomci( |
114 | 1x |
x = x2, n = n2, conf.level = conf.level, |
115 | 1x |
method = "wilsoncc" |
116 |
) |
|
117 | 1x |
l1 <- w1[2] |
118 | 1x |
u1 <- w1[3] |
119 | 1x |
l2 <- w2[2] |
120 | 1x |
u2 <- w2[3] |
121 | 1x |
ci_lwr <- max(-1, est - sqrt((p1_hat - l1)^2 + |
122 | 1x |
(u2 - p2_hat)^2)) |
123 | 1x |
ci_upr <- min(1, est + sqrt((u1 - p1_hat)^2 + (p2_hat - |
124 | 1x |
l2)^2)) |
125 |
}, |
|
126 | 18x |
mee = { |
127 | 1x |
.score <- function(p1, n1, p2, n2, dif) { |
128 | ! |
if (dif > 1) dif <- 1 |
129 | ! |
if (dif < -1) dif <- -1 |
130 | 24x |
diff <- p1 - p2 - dif |
131 | 24x |
if (abs(diff) == 0) { |
132 | ! |
res <- 0 |
133 |
} else { |
|
134 | 24x |
t <- n2 / n1 |
135 | 24x |
a <- 1 + t |
136 | 24x |
b <- -(1 + t + p1 + t * p2 + dif * (t + 2)) |
137 | 24x |
c <- dif * dif + dif * (2 * p1 + t + 1) + p1 + |
138 | 24x |
t * p2 |
139 | 24x |
d <- -p1 * dif * (1 + dif) |
140 | 24x |
v <- (b / a / 3)^3 - b * c / (6 * a * a) + d / a / 2 |
141 | 24x |
if (abs(v) < .Machine$double.eps) v <- 0 |
142 | 24x |
s <- sqrt((b / a / 3)^2 - c / a / 3) |
143 | 24x |
u <- ifelse(v > 0, 1, -1) * s |
144 | 24x |
w <- (3.141592654 + acos(v / u^3)) / 3 |
145 | 24x |
p1d <- 2 * u * cos(w) - b / a / 3 |
146 | 24x |
p2d <- p1d - dif |
147 | 24x |
n <- n1 + n2 |
148 | 24x |
res <- (p1d * (1 - p1d) / n1 + p2d * (1 - p2d) / n2) |
149 |
} |
|
150 | 24x |
return(sqrt(res)) |
151 |
} |
|
152 | 1x |
pval <- function(delta) { |
153 | 24x |
z <- (est - delta) / .score( |
154 | 24x |
p1_hat, n1, p2_hat, |
155 | 24x |
n2, delta |
156 |
) |
|
157 | 24x |
2 * min(stats::pnorm(z), 1 - stats::pnorm(z)) |
158 |
} |
|
159 | 1x |
ci_lwr <- max(-1, stats::uniroot(function(delta) { |
160 | 12x |
pval(delta) - |
161 | 12x |
alpha |
162 | 1x |
}, interval = c(-1 + 1e-06, est - 1e-06))$root) |
163 | 1x |
ci_upr <- min(1, stats::uniroot(function(delta) { |
164 | 12x |
pval(delta) - |
165 | 12x |
alpha |
166 | 1x |
}, interval = c(est + 1e-06, 1 - 1e-06))$root) |
167 |
}, |
|
168 | 18x |
blj = { |
169 | 1x |
p1_dash <- (x1 + 0.5) / (n1 + 1) |
170 | 1x |
p2_dash <- (x2 + 0.5) / (n2 + 1) |
171 | 1x |
vd <- p1_dash * (1 - p1_dash) / n1 + p2_dash * (1 - |
172 | 1x |
p2_dash) / n2 |
173 | 1x |
term2 <- kappa * sqrt(vd) |
174 | 1x |
est_dash <- p1_dash - p2_dash |
175 | 1x |
ci_lwr <- max(-1, est_dash - term2) |
176 | 1x |
ci_upr <- min(1, est_dash + term2) |
177 |
}, |
|
178 | 18x |
ha = { |
179 | 4x |
term2 <- 1 / (2 * min(n1, n2)) + kappa * sqrt(p1_hat * |
180 | 4x |
(1 - p1_hat) / (n1 - 1) + p2_hat * (1 - p2_hat) / (n2 - |
181 | 4x |
1)) |
182 | 4x |
ci_lwr <- max(-1, est - term2) |
183 | 4x |
ci_upr <- min(1, est + term2) |
184 |
}, |
|
185 | 18x |
mn = { |
186 | 1x |
.conf <- function(x1, n1, x2, n2, z, lower = FALSE) { |
187 | 2x |
p1 <- x1 / n1 |
188 | 2x |
p2 <- x2 / n2 |
189 | 2x |
p_hat <- p1 - p2 |
190 | 2x |
dp <- 1 + ifelse(lower, 1, -1) * p_hat |
191 | 2x |
i <- 1 |
192 | 2x |
while (i <= 50) { |
193 | 46x |
dp <- 0.5 * dp |
194 | 46x |
y <- p_hat + ifelse(lower, -1, 1) * dp |
195 | 46x |
score <- .score(p1, n1, p2, n2, y) |
196 | 46x |
if (score < z) { |
197 | 20x |
p_hat <- y |
198 |
} |
|
199 | 46x |
if ((dp < 1e-07) || (abs(z - score) < 1e-06)) { |
200 | 2x |
(break)() |
201 |
} else { |
|
202 | 44x |
i <- i + |
203 | 44x |
1 |
204 |
} |
|
205 |
} |
|
206 | 2x |
return(y) |
207 |
} |
|
208 | 1x |
.score <- function(p1, n1, p2, n2, dif) { |
209 | 46x |
diff <- p1 - p2 - dif |
210 | 46x |
if (abs(diff) == 0) { |
211 | ! |
res <- 0 |
212 |
} else { |
|
213 | 46x |
t <- n2 / n1 |
214 | 46x |
a <- 1 + t |
215 | 46x |
b <- -(1 + t + p1 + t * p2 + dif * (t + 2)) |
216 | 46x |
c <- dif * dif + dif * (2 * p1 + t + 1) + p1 + |
217 | 46x |
t * p2 |
218 | 46x |
d <- -p1 * dif * (1 + dif) |
219 | 46x |
v <- (b / a / 3)^3 - b * c / (6 * a * a) + d / a / 2 |
220 | 46x |
s <- sqrt((b / a / 3)^2 - c / a / 3) |
221 | 46x |
u <- ifelse(v > 0, 1, -1) * s |
222 | 46x |
w <- (3.141592654 + acos(v / u^3)) / 3 |
223 | 46x |
p1d <- 2 * u * cos(w) - b / a / 3 |
224 | 46x |
p2d <- p1d - dif |
225 | 46x |
n <- n1 + n2 |
226 | 46x |
var <- (p1d * (1 - p1d) / n1 + p2d * (1 - p2d) / n2) * |
227 | 46x |
n / (n - 1) |
228 | 46x |
res <- diff^2 / var |
229 |
} |
|
230 | 46x |
return(res) |
231 |
} |
|
232 | 1x |
z <- stats::qchisq(conf.level, 1) |
233 | 1x |
ci_lwr <- max(-1, .conf(x1, n1, x2, n2, z, TRUE)) |
234 | 1x |
ci_upr <- min(1, .conf(x1, n1, x2, n2, z, FALSE)) |
235 |
}, |
|
236 | 18x |
beal = { |
237 | ! |
a <- p1_hat + p2_hat |
238 | ! |
b <- p1_hat - p2_hat |
239 | ! |
u <- ((1 / n1) + (1 / n2)) / 4 |
240 | ! |
v <- ((1 / n1) - (1 / n2)) / 4 |
241 | ! |
V <- u * ((2 - a) * a - b^2) + 2 * v * (1 - a) * b # nolint |
242 | ! |
z <- stats::qchisq(p = 1 - alpha / 2, df = 1) |
243 | ! |
A <- sqrt(z * (V + z * u^2 * (2 - a) * a + z * v^2 * (1 - a)^2)) # nolint |
244 | ! |
B <- (b + z * v * (1 - a)) / (1 + z * u) # nolint |
245 | ! |
ci_lwr <- max(-1, B - A / (1 + z * u)) |
246 | ! |
ci_upr <- min(1, B + A / (1 + z * u)) |
247 |
}, |
|
248 | 18x |
hal = { |
249 | 1x |
psi <- (p1_hat + p2_hat) / 2 |
250 | 1x |
u <- (1 / n1 + 1 / n2) / 4 |
251 | 1x |
v <- (1 / n1 - 1 / n2) / 4 |
252 | 1x |
z <- kappa |
253 | 1x |
theta <- ((p1_hat - p2_hat) + z^2 * v * (1 - 2 * |
254 | 1x |
psi)) / (1 + z^2 * u) |
255 | 1x |
w <- z / (1 + z^2 * u) * sqrt(u * (4 * psi * (1 - psi) - |
256 | 1x |
(p1_hat - p2_hat)^2) + 2 * v * (1 - 2 * psi) * |
257 | 1x |
(p1_hat - p2_hat) + 4 * z^2 * u^2 * (1 - psi) * |
258 | 1x |
psi + z^2 * v^2 * (1 - 2 * psi)^2) |
259 | 1x |
c(theta + w, theta - w) |
260 | 1x |
ci_lwr <- max(-1, theta - w) |
261 | 1x |
ci_upr <- min(1, theta + w) |
262 |
}, |
|
263 | 18x |
jp = { |
264 | 1x |
psi <- 0.5 * ((x1 + 0.5) / (n1 + 1) + (x2 + 0.5) / (n2 + |
265 | 1x |
1)) |
266 | 1x |
u <- (1 / n1 + 1 / n2) / 4 |
267 | 1x |
v <- (1 / n1 - 1 / n2) / 4 |
268 | 1x |
z <- kappa |
269 | 1x |
theta <- ((p1_hat - p2_hat) + z^2 * v * (1 - 2 * |
270 | 1x |
psi)) / (1 + z^2 * u) |
271 | 1x |
w <- z / (1 + z^2 * u) * sqrt(u * (4 * psi * (1 - psi) - |
272 | 1x |
(p1_hat - p2_hat)^2) + 2 * v * (1 - 2 * psi) * |
273 | 1x |
(p1_hat - p2_hat) + 4 * z^2 * u^2 * (1 - psi) * |
274 | 1x |
psi + z^2 * v^2 * (1 - 2 * psi)^2) |
275 | 1x |
c(theta + w, theta - w) |
276 | 1x |
ci_lwr <- max(-1, theta - w) |
277 | 1x |
ci_upr <- min(1, theta + w) |
278 |
}, |
|
279 |
) |
|
280 | 18x |
ci <- c( |
281 | 18x |
est = est, lwr.ci = min(ci_lwr, ci_upr), |
282 | 18x |
upr.ci = max(ci_lwr, ci_upr) |
283 |
) |
|
284 | 18x |
if (sides == "left") { |
285 | ! |
ci[3] <- 1 |
286 | 18x |
} else if (sides == "right") { |
287 | ! |
ci[2] <- -1 |
288 |
} |
|
289 | 18x |
return(ci) |
290 |
} |
|
291 | 18x |
method <- match.arg(arg = method, several.ok = TRUE) |
292 | 18x |
sides <- match.arg(arg = sides, several.ok = TRUE) |
293 | 18x |
lst <- h_recycle( |
294 | 18x |
x1 = x1, n1 = n1, x2 = x2, n2 = n2, conf.level = conf.level, |
295 | 18x |
sides = sides, method = method |
296 |
) |
|
297 | 18x |
res <- t(sapply(1:attr(lst, "maxdim"), function(i) { |
298 | 18x |
iBinomDiffCI( |
299 | 18x |
x1 = lst$x1[i], |
300 | 18x |
n1 = lst$n1[i], x2 = lst$x2[i], n2 = lst$n2[i], conf.level = lst$conf.level[i], |
301 | 18x |
sides = lst$sides[i], method = lst$method[i] |
302 |
) |
|
303 |
})) |
|
304 | 18x |
lgn <- h_recycle(x1 = if (is.null(names(x1))) { |
305 | 18x |
paste("x1", seq_along(x1), sep = ".") |
306 |
} else { |
|
307 | ! |
names(x1) |
308 | 18x |
}, n1 = if (is.null(names(n1))) { |
309 | 18x |
paste("n1", seq_along(n1), sep = ".") |
310 |
} else { |
|
311 | ! |
names(n1) |
312 | 18x |
}, x2 = if (is.null(names(x2))) { |
313 | 18x |
paste("x2", seq_along(x2), sep = ".") |
314 |
} else { |
|
315 | ! |
names(x2) |
316 | 18x |
}, n2 = if (is.null(names(n2))) { |
317 | 18x |
paste("n2", seq_along(n2), sep = ".") |
318 |
} else { |
|
319 | ! |
names(n2) |
320 | 18x |
}, conf.level = conf.level, sides = sides, method = method) |
321 | 18x |
xn <- apply(as.data.frame(lgn[sapply(lgn, function(x) { |
322 | 126x |
length(unique(x)) != |
323 | 126x |
1 |
324 | 18x |
})]), 1, paste, collapse = ":") |
325 | 18x |
rownames(res) <- xn |
326 | 18x |
return(res) |
327 |
} |
|
328 | ||
329 |
#' @describeIn desctools_binom Compute confidence intervals for binomial proportions. |
|
330 |
#' |
|
331 |
#' @param x (`count`)\cr number of successes |
|
332 |
#' @param n (`count`)\cr number of trials |
|
333 |
#' @param conf.level (`proportion`)\cr confidence level, defaults to 0.95. |
|
334 |
#' @param sides (`character`)\cr side of the confidence interval to compute. Must be one of `"two-sided"` (default), |
|
335 |
#' `"left"`, or `"right"`. |
|
336 |
#' @param method (`character`)\cr method to use. Can be one out of: `"wald"`, `"wilson"`, `"wilsoncc"`, |
|
337 |
#' `"agresti-coull"`, `"jeffreys"`, `"modified wilson"`, `"modified jeffreys"`, `"clopper-pearson"`, `"arcsine"`, |
|
338 |
#' `"logit"`, `"witting"`, `"pratt"`, `"midp"`, `"lik"`, and `"blaker"`. |
|
339 |
#' |
|
340 |
#' @return A `matrix` with 3 columns containing: |
|
341 |
#' * `est`: estimate of proportion difference. |
|
342 |
#' * `lwr.ci`: lower end of the confidence interval. |
|
343 |
#' * `upr.ci`: upper end of the confidence interval. |
|
344 |
#' |
|
345 |
#' @keywords internal |
|
346 |
desctools_binomci <- function(x, |
|
347 |
n, |
|
348 |
conf.level = 0.95, # nolint |
|
349 |
sides = c("two.sided", "left", "right"), |
|
350 |
method = c( |
|
351 |
"wilson", "wald", "waldcc", "agresti-coull", |
|
352 |
"jeffreys", "modified wilson", "wilsoncc", "modified jeffreys", |
|
353 |
"clopper-pearson", "arcsine", "logit", "witting", "pratt", |
|
354 |
"midp", "lik", "blaker" |
|
355 |
), |
|
356 |
rand = 123, |
|
357 |
tol = 1e-05) { |
|
358 | 24x |
if (missing(method)) { |
359 | 1x |
method <- "wilson" |
360 |
} |
|
361 | 24x |
if (missing(sides)) { |
362 | 23x |
sides <- "two.sided" |
363 |
} |
|
364 | 24x |
iBinomCI <- function(x, n, conf.level = 0.95, sides = c( # nolint |
365 | 24x |
"two.sided", |
366 | 24x |
"left", "right" |
367 | 24x |
), method = c( |
368 | 24x |
"wilson", "wilsoncc", "wald", |
369 | 24x |
"waldcc", "agresti-coull", "jeffreys", "modified wilson", |
370 | 24x |
"modified jeffreys", "clopper-pearson", "arcsine", "logit", |
371 | 24x |
"witting", "pratt", "midp", "lik", "blaker" |
372 | 24x |
), rand = 123, |
373 | 24x |
tol = 1e-05) { |
374 | 24x |
if (length(x) != 1) { |
375 | ! |
stop("'x' has to be of length 1 (number of successes)") |
376 |
} |
|
377 | 24x |
if (length(n) != 1) { |
378 | ! |
stop("'n' has to be of length 1 (number of trials)") |
379 |
} |
|
380 | 24x |
if (length(conf.level) != 1) { |
381 | ! |
stop("'conf.level' has to be of length 1 (confidence level)") |
382 |
} |
|
383 | 24x |
if (conf.level < 0.5 || conf.level > 1) { |
384 | ! |
stop("'conf.level' has to be in [0.5, 1]") |
385 |
} |
|
386 | 24x |
sides <- match.arg(sides, choices = c( |
387 | 24x |
"two.sided", "left", |
388 | 24x |
"right" |
389 | 24x |
), several.ok = FALSE) |
390 | 24x |
if (sides != "two.sided") { |
391 | 1x |
conf.level <- 1 - 2 * (1 - conf.level) # nolint |
392 |
} |
|
393 | 24x |
alpha <- 1 - conf.level |
394 | 24x |
kappa <- stats::qnorm(1 - alpha / 2) |
395 | 24x |
p_hat <- x / n |
396 | 24x |
q_hat <- 1 - p_hat |
397 | 24x |
est <- p_hat |
398 | 24x |
switch(match.arg(arg = method, choices = c( |
399 | 24x |
"wilson", |
400 | 24x |
"wald", "waldcc", "wilsoncc", "agresti-coull", "jeffreys", |
401 | 24x |
"modified wilson", "modified jeffreys", "clopper-pearson", |
402 | 24x |
"arcsine", "logit", "witting", "pratt", "midp", "lik", |
403 | 24x |
"blaker" |
404 |
)), |
|
405 | 24x |
wald = { |
406 | 1x |
term2 <- kappa * sqrt(p_hat * q_hat) / sqrt(n) |
407 | 1x |
ci_lwr <- max(0, p_hat - term2) |
408 | 1x |
ci_upr <- min(1, p_hat + term2) |
409 |
}, |
|
410 | 24x |
waldcc = { |
411 | 1x |
term2 <- kappa * sqrt(p_hat * q_hat) / sqrt(n) |
412 | 1x |
term2 <- term2 + 1 / (2 * n) |
413 | 1x |
ci_lwr <- max(0, p_hat - term2) |
414 | 1x |
ci_upr <- min(1, p_hat + term2) |
415 |
}, |
|
416 | 24x |
wilson = { |
417 | 6x |
term1 <- (x + kappa^2 / 2) / (n + kappa^2) |
418 | 6x |
term2 <- kappa * sqrt(n) / (n + kappa^2) * sqrt(p_hat * |
419 | 6x |
q_hat + kappa^2 / (4 * n)) |
420 | 6x |
ci_lwr <- max(0, term1 - term2) |
421 | 6x |
ci_upr <- min(1, term1 + term2) |
422 |
}, |
|
423 | 24x |
wilsoncc = { |
424 | 3x |
lci <- (2 * x + kappa^2 - 1 - kappa * sqrt(kappa^2 - |
425 | 3x |
2 - 1 / n + 4 * p_hat * (n * q_hat + 1))) / (2 * |
426 | 3x |
(n + kappa^2)) |
427 | 3x |
uci <- (2 * x + kappa^2 + 1 + kappa * sqrt(kappa^2 + |
428 | 3x |
2 - 1 / n + 4 * p_hat * (n * q_hat - 1))) / (2 * |
429 | 3x |
(n + kappa^2)) |
430 | 3x |
ci_lwr <- max(0, ifelse(p_hat == 0, 0, lci)) |
431 | 3x |
ci_upr <- min(1, ifelse(p_hat == 1, 1, uci)) |
432 |
}, |
|
433 | 24x |
`agresti-coull` = { |
434 | 1x |
x_tilde <- x + kappa^2 / 2 |
435 | 1x |
n_tilde <- n + kappa^2 |
436 | 1x |
p_tilde <- x_tilde / n_tilde |
437 | 1x |
q_tilde <- 1 - p_tilde |
438 | 1x |
est <- p_tilde |
439 | 1x |
term2 <- kappa * sqrt(p_tilde * q_tilde) / sqrt(n_tilde) |
440 | 1x |
ci_lwr <- max(0, p_tilde - term2) |
441 | 1x |
ci_upr <- min(1, p_tilde + term2) |
442 |
}, |
|
443 | 24x |
jeffreys = { |
444 | 1x |
if (x == 0) { |
445 | ! |
ci_lwr <- 0 |
446 |
} else { |
|
447 | 1x |
ci_lwr <- stats::qbeta( |
448 | 1x |
alpha / 2, |
449 | 1x |
x + 0.5, n - x + 0.5 |
450 |
) |
|
451 |
} |
|
452 | 1x |
if (x == n) { |
453 | ! |
ci_upr <- 1 |
454 |
} else { |
|
455 | 1x |
ci_upr <- stats::qbeta(1 - |
456 | 1x |
alpha / 2, x + 0.5, n - x + 0.5) |
457 |
} |
|
458 |
}, |
|
459 | 24x |
`modified wilson` = { |
460 | 1x |
term1 <- (x + kappa^2 / 2) / (n + kappa^2) |
461 | 1x |
term2 <- kappa * sqrt(n) / (n + kappa^2) * sqrt(p_hat * |
462 | 1x |
q_hat + kappa^2 / (4 * n)) |
463 | 1x |
if ((n <= 50 & x %in% c(1, 2)) | (n >= 51 & x %in% |
464 | 1x |
c(1:3))) { |
465 | ! |
ci_lwr <- 0.5 * stats::qchisq(alpha, 2 * |
466 | ! |
x) / n |
467 |
} else { |
|
468 | 1x |
ci_lwr <- max(0, term1 - term2) |
469 |
} |
|
470 | 1x |
if ((n <= 50 & x %in% c(n - 1, n - 2)) | (n >= 51 & |
471 | 1x |
x %in% c(n - (1:3)))) { |
472 | ! |
ci_upr <- 1 - 0.5 * stats::qchisq( |
473 | ! |
alpha, |
474 | ! |
2 * (n - x) |
475 | ! |
) / n |
476 |
} else { |
|
477 | 1x |
ci_upr <- min(1, term1 + |
478 | 1x |
term2) |
479 |
} |
|
480 |
}, |
|
481 | 24x |
`modified jeffreys` = { |
482 | 1x |
if (x == n) { |
483 | ! |
ci_lwr <- (alpha / 2)^(1 / n) |
484 |
} else { |
|
485 | 1x |
if (x <= 1) { |
486 | ! |
ci_lwr <- 0 |
487 |
} else { |
|
488 | 1x |
ci_lwr <- stats::qbeta( |
489 | 1x |
alpha / 2, |
490 | 1x |
x + 0.5, n - x + 0.5 |
491 |
) |
|
492 |
} |
|
493 |
} |
|
494 | 1x |
if (x == 0) { |
495 | ! |
ci_upr <- 1 - (alpha / 2)^(1 / n) |
496 |
} else { |
|
497 | 1x |
if (x >= n - 1) { |
498 | ! |
ci_upr <- 1 |
499 |
} else { |
|
500 | 1x |
ci_upr <- stats::qbeta(1 - |
501 | 1x |
alpha / 2, x + 0.5, n - x + 0.5) |
502 |
} |
|
503 |
} |
|
504 |
}, |
|
505 | 24x |
`clopper-pearson` = { |
506 | 1x |
ci_lwr <- stats::qbeta(alpha / 2, x, n - x + 1) |
507 | 1x |
ci_upr <- stats::qbeta(1 - alpha / 2, x + 1, n - x) |
508 |
}, |
|
509 | 24x |
arcsine = { |
510 | 1x |
p_tilde <- (x + 0.375) / (n + 0.75) |
511 | 1x |
est <- p_tilde |
512 | 1x |
ci_lwr <- sin(asin(sqrt(p_tilde)) - 0.5 * kappa / sqrt(n))^2 |
513 | 1x |
ci_upr <- sin(asin(sqrt(p_tilde)) + 0.5 * kappa / sqrt(n))^2 |
514 |
}, |
|
515 | 24x |
logit = { |
516 | 1x |
lambda_hat <- log(x / (n - x)) |
517 | 1x |
V_hat <- n / (x * (n - x)) # nolint |
518 | 1x |
lambda_lower <- lambda_hat - kappa * sqrt(V_hat) |
519 | 1x |
lambda_upper <- lambda_hat + kappa * sqrt(V_hat) |
520 | 1x |
ci_lwr <- exp(lambda_lower) / (1 + exp(lambda_lower)) |
521 | 1x |
ci_upr <- exp(lambda_upper) / (1 + exp(lambda_upper)) |
522 |
}, |
|
523 | 24x |
witting = { |
524 | 1x |
set.seed(rand) |
525 | 1x |
x_tilde <- x + stats::runif(1, min = 0, max = 1) |
526 | 1x |
pbinom_abscont <- function(q, size, prob) { |
527 | 22x |
v <- trunc(q) |
528 | 22x |
term1 <- stats::pbinom(v - 1, size = size, prob = prob) |
529 | 22x |
term2 <- (q - v) * stats::dbinom(v, size = size, prob = prob) |
530 | 22x |
return(term1 + term2) |
531 |
} |
|
532 | 1x |
qbinom_abscont <- function(p, size, x) { |
533 | 2x |
fun <- function(prob, size, x, p) { |
534 | 22x |
pbinom_abscont(x, size, prob) - p |
535 |
} |
|
536 | 2x |
stats::uniroot(fun, |
537 | 2x |
interval = c(0, 1), size = size, |
538 | 2x |
x = x, p = p |
539 | 2x |
)$root |
540 |
} |
|
541 | 1x |
ci_lwr <- qbinom_abscont(1 - alpha, size = n, x = x_tilde) |
542 | 1x |
ci_upr <- qbinom_abscont(alpha, size = n, x = x_tilde) |
543 |
}, |
|
544 | 24x |
pratt = { |
545 | 1x |
if (x == 0) { |
546 | ! |
ci_lwr <- 0 |
547 | ! |
ci_upr <- 1 - alpha^(1 / n) |
548 | 1x |
} else if (x == 1) { |
549 | ! |
ci_lwr <- 1 - (1 - alpha / 2)^(1 / n) |
550 | ! |
ci_upr <- 1 - (alpha / 2)^(1 / n) |
551 | 1x |
} else if (x == (n - 1)) { |
552 | ! |
ci_lwr <- (alpha / 2)^(1 / n) |
553 | ! |
ci_upr <- (1 - alpha / 2)^(1 / n) |
554 | 1x |
} else if (x == n) { |
555 | ! |
ci_lwr <- alpha^(1 / n) |
556 | ! |
ci_upr <- 1 |
557 |
} else { |
|
558 | 1x |
z <- stats::qnorm(1 - alpha / 2) |
559 | 1x |
A <- ((x + 1) / (n - x))^2 # nolint |
560 | 1x |
B <- 81 * (x + 1) * (n - x) - 9 * n - 8 # nolint |
561 | 1x |
C <- (0 - 3) * z * sqrt(9 * (x + 1) * (n - x) * (9 * n + 5 - z^2) + n + 1) # nolint |
562 | 1x |
D <- 81 * (x + 1)^2 - 9 * (x + 1) * (2 + z^2) + 1 # nolint |
563 | 1x |
E <- 1 + A * ((B + C) / D)^3 # nolint |
564 | 1x |
ci_upr <- 1 / E |
565 | 1x |
A <- (x / (n - x - 1))^2 # nolint |
566 | 1x |
B <- 81 * x * (n - x - 1) - 9 * n - 8 # nolint |
567 | 1x |
C <- 3 * z * sqrt(9 * x * (n - x - 1) * (9 * n + 5 - z^2) + n + 1) # nolint |
568 | 1x |
D <- 81 * x^2 - 9 * x * (2 + z^2) + 1 # nolint |
569 | 1x |
E <- 1 + A * ((B + C) / D)^3 # nolint |
570 | 1x |
ci_lwr <- 1 / E |
571 |
} |
|
572 |
}, |
|
573 | 24x |
midp = { |
574 | 1x |
f_low <- function(pi, x, n) { |
575 | 12x |
1 / 2 * stats::dbinom(x, size = n, prob = pi) + stats::pbinom(x, |
576 | 12x |
size = n, prob = pi, lower.tail = FALSE |
577 |
) - |
|
578 | 12x |
(1 - conf.level) / 2 |
579 |
} |
|
580 | 1x |
f_up <- function(pi, x, n) { |
581 | 12x |
1 / 2 * stats::dbinom(x, size = n, prob = pi) + stats::pbinom(x - |
582 | 12x |
1, size = n, prob = pi) - (1 - conf.level) / 2 |
583 |
} |
|
584 | 1x |
ci_lwr <- 0 |
585 | 1x |
ci_upr <- 1 |
586 | 1x |
if (x != 0) { |
587 | 1x |
ci_lwr <- stats::uniroot(f_low, |
588 | 1x |
interval = c(0, p_hat), |
589 | 1x |
x = x, n = n |
590 | 1x |
)$root |
591 |
} |
|
592 | 1x |
if (x != n) { |
593 | 1x |
ci_upr <- stats::uniroot(f_up, interval = c( |
594 | 1x |
p_hat, |
595 | 1x |
1 |
596 | 1x |
), x = x, n = n)$root |
597 |
} |
|
598 |
}, |
|
599 | 24x |
lik = { |
600 | 2x |
ci_lwr <- 0 |
601 | 2x |
ci_upr <- 1 |
602 | 2x |
z <- stats::qnorm(1 - alpha * 0.5) |
603 | 2x |
tol <- .Machine$double.eps^0.5 |
604 | 2x |
BinDev <- function(y, x, mu, wt, bound = 0, tol = .Machine$double.eps^0.5, # nolint |
605 |
...) { |
|
606 | 40x |
ll_y <- ifelse(y %in% c(0, 1), 0, stats::dbinom(x, wt, |
607 | 40x |
y, |
608 | 40x |