diff --git a/.Rbuildignore b/.Rbuildignore index f16f2b48..95533ebb 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -22,3 +22,4 @@ init.sh workflows.md images ^data-raw$ +^\.gitleaksignore$ diff --git a/.github/ISSUE_TEMPLATE/bug.yml b/.github/ISSUE_TEMPLATE/bug.yml new file mode 100644 index 00000000..9808a314 --- /dev/null +++ b/.github/ISSUE_TEMPLATE/bug.yml @@ -0,0 +1,50 @@ +--- +name: 🐞 Bug Report +description: File a bug report +title: "[Bug]: " +labels: ["bug"] +assignees: + - danielinteractive +body: + - type: markdown + attributes: + value: | + Thanks for taking the time to fill out this bug report! + - type: textarea + id: what-happened + attributes: + label: What happened? + description: Also tell us, what did you expect to happen? + placeholder: Tell us what you see! + value: "A bug happened!" + validations: + required: true + - type: textarea + id: logs + attributes: + label: Relevant log output + description: Please copy and paste any relevant log output. This will be automatically formatted into code, so no need for backticks. + - type: checkboxes + id: code-of-conduct + attributes: + label: Code of Conduct + description: By submitting this issue, you agree to follow our [Code of Conduct.](https://github.com/insightsengineering/.github/blob/main/CODE_OF_CONDUCT.md) + options: + - label: I agree to follow this project's Code of Conduct. + required: true + - type: checkboxes + id: contributor-guidelines + attributes: + label: Contribution Guidelines + description: By submitting this issue, you agree to follow our [Contribution Guidelines.](https://github.com/insightsengineering/.github/blob/main/CONTRIBUTING.md) + options: + - label: I agree to follow this project's Contribution Guidelines. + required: true + - type: checkboxes + id: security-policy + attributes: + label: Security Policy + description: By submitting this issue, you agree to follow our [Security Policy.](https://github.com/insightsengineering/.github/security/policy) + options: + - label: I agree to follow this project's Security Policy. + required: true diff --git a/.github/ISSUE_TEMPLATE/config.yml b/.github/ISSUE_TEMPLATE/config.yml new file mode 100644 index 00000000..8612ecdb --- /dev/null +++ b/.github/ISSUE_TEMPLATE/config.yml @@ -0,0 +1,10 @@ +--- +blank_issues_enabled: true + +contact_links: + - name: We are hiring! + url: https://careers.gene.com/ + about: Genentech and Roche are hiring! + - name: openstatsware + url: https://openstatsware.org/ + about: Software Engineering in Biostatistics diff --git a/.github/ISSUE_TEMPLATE/feature.yml b/.github/ISSUE_TEMPLATE/feature.yml new file mode 100644 index 00000000..934e269a --- /dev/null +++ b/.github/ISSUE_TEMPLATE/feature.yml @@ -0,0 +1,37 @@ +--- +name: ✹ Feature Request +description: Request or propose a new feature +title: "[Feature Request]: <title>" +labels: ["enhancement"] +assignees: + - danielinteractive +body: + - type: textarea + attributes: + label: Feature description + validations: + required: true + - type: checkboxes + id: code-of-conduct + attributes: + label: Code of Conduct + description: By submitting this issue, you agree to follow our [Code of Conduct.](https://github.com/insightsengineering/.github/blob/main/CODE_OF_CONDUCT.md) + options: + - label: I agree to follow this project's Code of Conduct. + required: true + - type: checkboxes + id: contributor-guidelines + attributes: + label: Contribution Guidelines + description: By submitting this issue, you agree to follow our [Contribution Guidelines.](https://github.com/insightsengineering/.github/blob/main/CONTRIBUTING.md) + options: + - label: I agree to follow this project's Contribution Guidelines. + required: true + - type: checkboxes + id: security-policy + attributes: + label: Security Policy + description: By submitting this issue, you agree to follow our [Security Policy.](https://github.com/insightsengineering/.github/security/policy) + options: + - label: I agree to follow this project's Security Policy. + required: true diff --git a/.github/ISSUE_TEMPLATE/question.yml b/.github/ISSUE_TEMPLATE/question.yml new file mode 100644 index 00000000..4afdaeff --- /dev/null +++ b/.github/ISSUE_TEMPLATE/question.yml @@ -0,0 +1,37 @@ +--- +name: ❓ Question +description: Question about usage or documentation +title: "[Question]: <title>" +labels: ["question"] +assignees: + - danielinteractive +body: + - type: textarea + attributes: + label: What is your question? + validations: + required: true + - type: checkboxes + id: code-of-conduct + attributes: + label: Code of Conduct + description: By submitting this issue, you agree to follow our [Code of Conduct.](https://github.com/insightsengineering/.github/blob/main/CODE_OF_CONDUCT.md) + options: + - label: I agree to follow this project's Code of Conduct. + required: true + - type: checkboxes + id: contributor-guidelines + attributes: + label: Contribution Guidelines + description: By submitting this issue, you agree to follow our [Contribution Guidelines.](https://github.com/insightsengineering/.github/blob/main/CONTRIBUTING.md) + options: + - label: I agree to follow this project's Contribution Guidelines. + required: true + - type: checkboxes + id: security-policy + attributes: + label: Security Policy + description: By submitting this issue, you agree to follow our [Security Policy.](https://github.com/insightsengineering/.github/security/policy) + options: + - label: I agree to follow this project's Security Policy. + required: true diff --git a/.gitleaksignore b/.gitleaksignore new file mode 100644 index 00000000..4c7ddefc --- /dev/null +++ b/.gitleaksignore @@ -0,0 +1,2 @@ +R/RcppExports.R:generic-api-key:2 +src/RcppExports.cpp:generic-api-key:2 diff --git a/DESCRIPTION b/DESCRIPTION index 8ffc4939..15ee3d4a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,10 +1,13 @@ Type: Package Package: shrinkforest -Title: Shrinkage based Forest Plots +Title: Shrinkage Based Forest Plots Version: 0.0.0.9033 -Date: 2023-02-28 +Date: 2024-03-07 Authors@R: - person("Mar", "Vazquez Rabunal", , "mar.vazquez_rabunal@roche.com", role = c("aut", "cre")) + c( + person("Mar", "Vazquez Rabunal", , "mar.vazquez_rabunal@roche.com", role = c("aut", "cre")), + person("Daniel", "SabanĂ©s BovĂ©", , "daniel.sabanes_bove@roche.com", role = "aut") + ) Description: Subgroup analyses are routinely performed in clinical trial analyses. From a methodological perspective, two key issues of subgroup analyses are multiplicity (even if only predefined subgroups @@ -49,4 +52,4 @@ Language: en-US LazyData: true LazyDataCompression: xz Roxygen: list(markdown = TRUE) -RoxygenNote: 7.2.3 +RoxygenNote: 7.3.1 diff --git a/R/generate_stacked_data.R b/R/generate_stacked_data.R index e8ebe55d..0789fcd8 100644 --- a/R/generate_stacked_data.R +++ b/R/generate_stacked_data.R @@ -37,18 +37,24 @@ generate_stacked_data <- function(base_model, subgroup_model, data, } tmp <- lapply(data[, subgroup_vars], function(x) levels(factor(x))) subgroup_names <- paste(rep(names(tmp), lapply(tmp, length)), - unlist(tmp), sep = ".") + unlist(tmp), + sep = "." + ) data <- data[, c(arm_var, resp_var, status_var, subgroup_vars)] data[, subgroup_vars] <- lapply(data[, subgroup_vars], as.character) - d <- tidyr::gather(data, "subgroup_var", "subgroup_value", - -arm_var, -resp_var, -status_var) + d <- tidyr::gather( + data, "subgroup_var", "subgroup_value", + -arm_var, -resp_var, -status_var + ) subgroup <- paste(d$subgroup_var, d$subgroup_value, sep = ".") subgroup <- factor(subgroup, levels = subgroup_names) d <- cbind(d, subgroup) d <- dplyr::arrange(d, subgroup) d <- if (resptype == "survival") { - dplyr::rename_at(d, dplyr::vars(c(arm_var, resp_var, status_var)), - ~ c("arm", "time", "status")) + dplyr::rename_at( + d, dplyr::vars(c(arm_var, resp_var, status_var)), + ~ c("arm", "time", "status") + ) } else if (resptype == "binary") { dplyr::rename_at(d, dplyr::vars(c(arm_var, resp_var)), ~ c("arm", "y")) } diff --git a/R/horseshoe.R b/R/horseshoe.R index b96164ae..7ff3c312 100644 --- a/R/horseshoe.R +++ b/R/horseshoe.R @@ -29,7 +29,8 @@ #' @examples #' horseshoe("ev_pfs", "arm", c("x_1", "x_2"), c("x_1", "x_2", "x_3"), #' example_data, "binary", -#' chains = 1, seed = 0, control = list(adapt_delta = 0.95) +#' chains = 1, seed = 0, control = list(adapt_delta = 0.95), +#' iter = 50 # In practice, you need to omit this or set it much higher! #' ) horseshoe <- function(resp, trt, subgr, covars, data, resptype = c("survival", "binary"), status = NULL, @@ -75,9 +76,12 @@ horseshoe <- function(resp, trt, subgr, covars, data, design_main <- prep_data$design_main[, -1] design_matrix <- cbind(design_main, prep_data$design_ia) form <- stats::as.formula(paste(resp, " ~ a + b")) - form_a <- stats::as.formula(paste("a ~ 1 +", paste0(colnames(design_main), - collapse = " + " - ))) + form_a <- stats::as.formula( + paste( + "a ~ 1 +", + paste0(colnames(design_main), collapse = " + ") + ) + ) y <- as.data.frame(data[[resp]]) colnames(y) <- resp data_model <- cbind(design_matrix, y) @@ -90,11 +94,8 @@ horseshoe <- function(resp, trt, subgr, covars, data, brms::lf(form_b), data = data_model, family = family, - brms::prior(normal(0, 5), class = "b", nlpar = "a") + - brms::prior(horseshoe(1), - class = "b", - nlpar = "b" - ), + brms::prior_string("normal(0, 5)", class = "b", nlpar = "a") + + brms::prior_string("horseshoe(1)", class = "b", nlpar = "b"), ... ) } diff --git a/R/lOR_estimation.R b/R/lOR_estimation.R index 0c582db5..af9d7da3 100644 --- a/R/lOR_estimation.R +++ b/R/lOR_estimation.R @@ -18,18 +18,22 @@ lor_estimation <- function(x_subg, dummy_subg, est_coef) { assert_matrix(dummy_subg) assert_matrix(est_coef) x_arm <- list() - for (i in 0:1){ - x_arm[[i + 1]] <- cbind(rep(1, nrow(x_subg)), rep(i, nrow(x_subg)), - x_subg[, 2:(ncol(x_subg) - ncol(dummy_subg))], i * dummy_subg) + for (i in 0:1) { + x_arm[[i + 1]] <- cbind( + rep(1, nrow(x_subg)), rep(i, nrow(x_subg)), + x_subg[, 2:(ncol(x_subg) - ncol(dummy_subg))], i * dummy_subg + ) } prob <- function(x, est_coef) { - k <- as.matrix(x) %*% as.matrix(est_coef) + k <- as.matrix(x) %*% as.matrix(est_coef) p <- exp(k) / (1 + exp(k)) y <- apply(p, 2, mean) y } y_arm <- lapply(x_arm, prob, est_coef = est_coef) - phi <- as.numeric(log(y_arm[[2]] / (1 - y_arm[[2]])) - - log(y_arm[[1]] / (1 - y_arm[[1]]))) + phi <- as.numeric( + log(y_arm[[2]] / (1 - y_arm[[2]])) - + log(y_arm[[1]] / (1 - y_arm[[1]])) + ) phi } diff --git a/R/plot.R b/R/plot.R index 430dc858..a80be7b4 100644 --- a/R/plot.R +++ b/R/plot.R @@ -16,12 +16,24 @@ plot.summary.naive <- function(x, ...) { data$subgroup <- as.factor(data$subgroup) forestplot <- ggplot( data = data, - aes(x = trt.estimate, y = model, xmin = trt.low, xmax = trt.high) + aes( + x = .data$trt.estimate, + y = .data$model, + xmin = .data$trt.low, + xmax = .data$trt.high + ) ) + ggtitle("Forest plot Naive") + - geom_pointrange(aes(col = model)) + + geom_pointrange(aes(col = .data$model)) + ylab("Subgroup") + - geom_errorbar(aes(xmin = trt.low, xmax = trt.high, col = model), width = 0.5, cex = 1) + + geom_errorbar( + aes( + xmin = .data$trt.low, + xmax = .data$trt.high, + col = .data$model + ), + width = 0.5, cex = 1 + ) + facet_wrap(~ forcats::fct_inorder(subgroup), strip.position = "left", nrow = nrow(data), scales = "free_y" @@ -70,10 +82,10 @@ plot.summary.elastic_net <- function(x, ...) { data$subgroup <- as.factor(data$subgroup) forestplot <- ggplot( data = data, - aes(x = trt.estimate, y = model) + aes(x = .data$trt.estimate, y = .data$model) ) + ggtitle("Forest plot Elastic Net") + - geom_point(aes(col = model)) + + geom_point(aes(col = .data$model)) + ylab("Subgroup") + facet_wrap(~ forcats::fct_inorder(subgroup), strip.position = "left", @@ -117,12 +129,25 @@ plot.summary.horseshoe <- function(x, ...) { data$subgroup <- as.factor(data$subgroup) forestplot <- ggplot( data = data, - aes(x = trt.estimate, y = model, xmin = trt.low, xmax = trt.high) + aes( + x = .data$trt.estimate, + y = .data$model, + xmin = .data$trt.low, + xmax = .data$trt.high + ) ) + ggtitle("Forest plot Horseshoe") + - geom_pointrange(aes(col = model)) + + geom_pointrange(aes(col = .data$model)) + ylab("Subgroup") + - geom_errorbar(aes(xmin = trt.low, xmax = trt.high, col = model), width = 0.5, cex = 1) + + geom_errorbar( + aes( + xmin = .data$trt.low, + xmax = .data$trt.high, + col = .data$model + ), + width = 0.5, + cex = 1 + ) + facet_wrap(~ forcats::fct_inorder(subgroup), strip.position = "left", nrow = nrow(data), scales = "free_y" @@ -170,16 +195,19 @@ plot.compare.data <- function(x, ...) { forestplot <- ggplot( data = data, aes( - x = trt.estimate, y = forcats::fct_rev(forcats::fct_inorder(model)), xmin = trt.low, - xmax = trt.high + x = .data$trt.estimate, + y = forcats::fct_rev(forcats::fct_inorder(.data$model)), + xmin = .data$trt.low, + xmax = .data$trt.high ) ) + ggtitle("Forest plot") + - geom_pointrange(aes(col = forcats::fct_inorder(model))) + + geom_pointrange(aes(col = forcats::fct_inorder(.data$model))) + ylab("Subgroup") + geom_errorbar(aes( - xmin = trt.low, xmax = trt.high, - col = forcats::fct_inorder(model) + xmin = .data$trt.low, + xmax = .data$trt.high, + col = forcats::fct_inorder(.data$model) ), width = 0.5, cex = 1) + facet_wrap(~ forcats::fct_inorder(subgroup), strip.position = "left", diff --git a/R/trt_horseshoe.R b/R/trt_horseshoe.R index 1b9242b1..a35fb57e 100644 --- a/R/trt_horseshoe.R +++ b/R/trt_horseshoe.R @@ -19,7 +19,7 @@ #' @export #' #' @examples -#' trt_horseshoe(horseshoe_fit_surv) +#' trt_horseshoe(horseshoe_fit_surv, m = 5) trt_horseshoe <- function(object, gamma = 1, l = NULL, m = 50) { assert_class(object, c("shrinkforest", "horseshoe")) assert_int(m) diff --git a/data-raw/elastic_net_fit_surv.R b/data-raw/elastic_net_fit_surv.R index fb6e241d..1651bd74 100644 --- a/data-raw/elastic_net_fit_surv.R +++ b/data-raw/elastic_net_fit_surv.R @@ -1,7 +1,13 @@ -elastic_net_fit_surv <- elastic_net("tt_pfs", "arm", c("x_1", "x_2", "x_3", - "x_4", "x_5", "x_6", "x_7", - "x_8", "x_9", "x_10"), - c("x_1", "x_2", "x_3", "x_4", "x_5", "x_6", - "x_7", "x_8", "x_9", "x_10"), - example_data, "survival", 1, "ev_pfs") +elastic_net_fit_surv <- elastic_net( + "tt_pfs", "arm", c( + "x_1", "x_2", "x_3", + "x_4", "x_5", "x_6", "x_7", + "x_8", "x_9", "x_10" + ), + c( + "x_1", "x_2", "x_3", "x_4", "x_5", "x_6", + "x_7", "x_8", "x_9", "x_10" + ), + example_data, "survival", 1, "ev_pfs" +) usethis::use_data(elastic_net_fit_surv, compress = "xz", overwrite = TRUE) diff --git a/data-raw/est_coef_bin1.R b/data-raw/est_coef_bin1.R index 0c3ae0ca..c9759d9d 100644 --- a/data-raw/est_coef_bin1.R +++ b/data-raw/est_coef_bin1.R @@ -1,3 +1,4 @@ est_coef_bin1 <- as.matrix(coef(elastic_net_fit_bin$fit, - s = elastic_net_fit_bin$fit$lambda.min)) + s = elastic_net_fit_bin$fit$lambda.min +)) usethis::use_data(est_coef_bin1, overwrite = TRUE) diff --git a/data-raw/example_data.R b/data-raw/example_data.R index 0da645ed..429de607 100644 --- a/data-raw/example_data.R +++ b/data-raw/example_data.R @@ -9,48 +9,68 @@ simul_covariates <- function(n, p_catvar = 10, add_contvars = FALSE) { # block diagonal covariance matrix for underlying multivariate normal data # create covariate matrix in blocks of 10 sigma <- matrix(0, nrow = 10, ncol = 10) - sigma[1:5, 1:5] <- 0 # first 5 covariates uncorrelated with everything - sigma[6:8, 6:8] <- 0.25 # cov 6-8 with "moderate" correlation - sigma[9:10, 9:10] <- 0.5 # cov 9-10 with "high" correlation + sigma[1:5, 1:5] <- 0 # first 5 covariates uncorrelated with everything + sigma[6:8, 6:8] <- 0.25 # cov 6-8 with "moderate" correlation + sigma[9:10, 9:10] <- 0.5 # cov 9-10 with "high" correlation diag(sigma) <- 1 # variance 1 no_10_blocks <- ceiling(p_catvar / 10) x <- NULL z <- NULL - for (j in 1:no_10_blocks){ + for (j in 1:no_10_blocks) { # continuous version z_j <- data.frame(mvrnorm(n, mu = rep(0, 10), Sigma = sigma)) colnames(z_j) <- paste("z", (j - 1) * 10 + 1:10, sep = "_") - if (j == 1) {z <- z_j} else { z <- cbind(z, z_j)} #nolint + if (j == 1) { + z <- z_j + } else { + z <- cbind(z, z_j) + } # nolint # categorized version - x_j <- data.frame(v1 = cut(z_j[, 1], c(-Inf, qnorm(0.5), Inf), - labels = c("a", "b")), - v2 = cut(z_j[, 2], c(-Inf, qnorm(0.4), Inf), - labels = c("a", "b")), - v3 = cut(z_j[, 3], c(-Inf, qnorm(0.2), Inf), - labels = c("a", "b")), - v4 = cut(z_j[, 4], c(-Inf, qnorm(c(0.3, 0.6)), Inf), - labels = c("a", "b", "c")), - v5 = cut(z_j[, 5], c(-Inf, qnorm(c(0.15, 0.3, 0.6)), Inf), - labels = c("a", "b", "c", "d")), - v6 = cut(z_j[, 6], c(-Inf, qnorm(0.4), Inf), - labels = c("a", "b")), - v7 = cut(z_j[, 7], c(-Inf, qnorm(0.4), Inf), - labels = c("a", "b")), - v8 = cut(z_j[, 8], c(-Inf, qnorm(c(0.2, 0.5)), Inf), - labels = c("a", "b", "c")), - v9 = cut(z_j[, 9], c(-Inf, qnorm(0.2), Inf), - labels = c("a", "b")), - v10 = cut(z_j[, 10], c(-Inf, qnorm(c(0.2, 0.5)), Inf), - labels = c("a", "b", "c"))) + x_j <- data.frame( + v1 = cut(z_j[, 1], c(-Inf, qnorm(0.5), Inf), + labels = c("a", "b") + ), + v2 = cut(z_j[, 2], c(-Inf, qnorm(0.4), Inf), + labels = c("a", "b") + ), + v3 = cut(z_j[, 3], c(-Inf, qnorm(0.2), Inf), + labels = c("a", "b") + ), + v4 = cut(z_j[, 4], c(-Inf, qnorm(c(0.3, 0.6)), Inf), + labels = c("a", "b", "c") + ), + v5 = cut(z_j[, 5], c(-Inf, qnorm(c(0.15, 0.3, 0.6)), Inf), + labels = c("a", "b", "c", "d") + ), + v6 = cut(z_j[, 6], c(-Inf, qnorm(0.4), Inf), + labels = c("a", "b") + ), + v7 = cut(z_j[, 7], c(-Inf, qnorm(0.4), Inf), + labels = c("a", "b") + ), + v8 = cut(z_j[, 8], c(-Inf, qnorm(c(0.2, 0.5)), Inf), + labels = c("a", "b", "c") + ), + v9 = cut(z_j[, 9], c(-Inf, qnorm(0.2), Inf), + labels = c("a", "b") + ), + v10 = cut(z_j[, 10], c(-Inf, qnorm(c(0.2, 0.5)), Inf), + labels = c("a", "b", "c") + ) + ) colnames(x_j) <- paste("x", (j - 1) * 10 + 1:10, sep = "_") - if (j == 1) {x <- x_j} else { x <- cbind(x, x_j)} #nolint + if (j == 1) { + x <- x_j + } else { + x <- cbind(x, x_j) + } # nolint } x <- cbind(arm = sample(rep(c(0, 1), c(n %/% 2, n - n %/% 2))), x[, 1:p_catvar]) if (add_contvars) x <- cbind(x, z[, 1:p_catvar]) x } -simul_pfs <- function(lp_aft, sigma_aft, recr_duration, rate_cens, n_events){ #nolint +simul_pfs <- function(lp_aft, sigma_aft, recr_duration, rate_cens, n_events) { # nolint n <- length(lp_aft) # Uncensored event time log_tt_pfs <- c(lp_aft + sigma_aft * log(rexp(n, rate = 1))) @@ -59,34 +79,39 @@ simul_pfs <- function(lp_aft, sigma_aft, recr_duration, rate_cens, n_events){ #n tt_pfs_cens1 <- rexp(n, rate = rate_cens) tt_pfs_cens1 <- pmin(tt_pfs_uncens, tt_pfs_cens1) ev_pfs_cens1 <- ifelse(tt_pfs_uncens <= tt_pfs_cens1, 1, 0) - if (sum(ev_pfs_cens1) < n_events) stop(paste("Impossible to reach", n_events, - "events with", n, "patients, + if (sum(ev_pfs_cens1) < n_events) { + stop(paste( + "Impossible to reach", n_events, + "events with", n, "patients, a censoring rate of", rate_cens, - "and the specified linear predictor.")) + "and the specified linear predictor." + )) + } # censoring 2: due to staggerred recruitment and recruiting only until target_ev events have been observed rec_time <- runif(n, min = 0, max = recr_duration) tt_pfs_cens1_calendar <- rec_time + tt_pfs_cens1 study_stop_time <- sort(tt_pfs_cens1_calendar[ev_pfs_cens1 == 1])[n_events] - if (study_stop_time < max(rec_time)) + if (study_stop_time < max(rec_time)) { warning("Target number of events reached before all subjects were enrolled.") + } tt_pfs <- pmax(0, pmin(tt_pfs_cens1_calendar, study_stop_time) - rec_time) ev_pfs <- ifelse(tt_pfs_cens1_calendar <= study_stop_time, ev_pfs_cens1, 0) data.frame(tt_pfs = tt_pfs, ev_pfs = ev_pfs) } -quicksimul <- function(n, coef, sigma_aft, recr_duration, rate_cens, n_events){ #nolint +quicksimul <- function(n, coef, sigma_aft, recr_duration, rate_cens, n_events) { # nolint # Quickly simulate actual data combining functions covariates - #and simul_pfs (assuming 10 covariates) + # and simul_pfs (assuming 10 covariates) # Regression coefficients are for an AFT with over-parametrized dummy - #coding for arm-subgroup interactions (see creation of design matrix below) + # coding for arm-subgroup interactions (see creation of design matrix below) covariates <- simul_covariates(n = n, p_catvar = 10, add_contvars = FALSE) #- create design matrix with over-parametrized - #dummy coding for arm-subgroup interactions + # dummy coding for arm-subgroup interactions subgroup_model <- ~ x_1 + x_2 + x_3 + x_4 + x_5 + x_6 + x_7 + x_8 + x_9 + x_10 design_main <- model.matrix(update(subgroup_model, ~ arm + .), data = covariates) subgroup_vars <- all.vars(subgroup_model) design_ia <- NULL - for (j in subgroup_vars){ + for (j in subgroup_vars) { ia_j <- model.matrix(as.formula(paste("~", j, "-1")), data = covariates) * covariates$arm design_ia <- cbind(design_ia, ia_j) } @@ -95,30 +120,34 @@ quicksimul <- function(n, coef, sigma_aft, recr_duration, rate_cens, n_events){ design_matrix <- cbind(design_main, design_ia) #- get linear predictor for AFT and simulate corresponding outcome lp_aft <- design_matrix %*% coef # linear predictor - outcome <- simul_pfs(lp_aft = lp_aft, sigma_aft = sigma_aft, - recr_duration = recr_duration, rate_cens = rate_cens, - n_events = n_events) + outcome <- simul_pfs( + lp_aft = lp_aft, sigma_aft = sigma_aft, + recr_duration = recr_duration, rate_cens = rate_cens, + n_events = n_events + ) d <- cbind(id = 1:n, covariates, outcome) d } -#Set seed +# Set seed set.seed(0) # sigma_aft sigma_aft <- 0.85 # Dummy vector for AFT regression coefficient (with over-parametrized dummy -#coding for arm-covariate interactions) +# coding for arm-covariate interactions) coef_raw <- rep(0, 42) -names(coef_raw) <- c("(Intercept)", "arm", "x_1b", "x_2b", "x_3b", "x_4b", "x_4c", - "x_5b", "x_5c", "x_5d", "x_6b", "x_7b", "x_8b", "x_8c", - "x_9b", "x_10b", "x_10c", "x_1a_arm", "x_1b_arm", "x_2a_arm", - "x_2b_arm", "x_3a_arm", "x_3b_arm", "x_4a_arm", "x_4b_arm", - "x_4c_arm", "x_5a_arm", "x_5b_arm", "x_5c_arm", "x_5d_arm", - "x_6a_arm", "x_6b_arm", "x_7a_arm", "x_7b_arm", "x_8a_arm", - "x_8b_arm", "x_8c_arm", "x_9a_arm", "x_9b_arm", "x_10a_arm", - "x_10b_arm", "x_10c_arm") +names(coef_raw) <- c( + "(Intercept)", "arm", "x_1b", "x_2b", "x_3b", "x_4b", "x_4c", + "x_5b", "x_5c", "x_5d", "x_6b", "x_7b", "x_8b", "x_8c", + "x_9b", "x_10b", "x_10c", "x_1a_arm", "x_1b_arm", "x_2a_arm", + "x_2b_arm", "x_3a_arm", "x_3b_arm", "x_4a_arm", "x_4b_arm", + "x_4c_arm", "x_5a_arm", "x_5b_arm", "x_5c_arm", "x_5d_arm", + "x_6a_arm", "x_6b_arm", "x_7a_arm", "x_7b_arm", "x_8a_arm", + "x_8b_arm", "x_8c_arm", "x_9a_arm", "x_9b_arm", "x_10a_arm", + "x_10b_arm", "x_10c_arm" +) # Set intercept and prognostic factors (same for all scenarios) @@ -132,7 +161,9 @@ coef_raw["x_6b"] <- -log(1.5) * sigma_aft coef_1 <- coef_raw coef_1["arm"] <- -log(0.66) * sigma_aft -example_data <- quicksimul(n = 1000, coef = coef_1, sigma_aft = sigma_aft, - recr_duration = 3, rate_cens = 0.02, n_events = 247) +example_data <- quicksimul( + n = 1000, coef = coef_1, sigma_aft = sigma_aft, + recr_duration = 3, rate_cens = 0.02, n_events = 247 +) example_data$arm <- as.factor(example_data$arm) usethis::use_data(example_data) diff --git a/data-raw/naive_fit_bin.R b/data-raw/naive_fit_bin.R index 2d040160..b68984b4 100644 --- a/data-raw/naive_fit_bin.R +++ b/data-raw/naive_fit_bin.R @@ -1,4 +1,8 @@ -naive_fit_bin <- naive("ev_pfs", "arm", c("x_1", "x_2", "x_3", "x_4", "x_5", - "x_6", "x_7", "x_8", "x_9", "x_10"), - example_data, "binary") +naive_fit_bin <- naive( + "ev_pfs", "arm", c( + "x_1", "x_2", "x_3", "x_4", "x_5", + "x_6", "x_7", "x_8", "x_9", "x_10" + ), + example_data, "binary" +) usethis::use_data(naive_fit_bin, compress = "xz", overwrite = TRUE) diff --git a/data-raw/naive_fit_surv.R b/data-raw/naive_fit_surv.R index babf0154..9d0c3ca4 100644 --- a/data-raw/naive_fit_surv.R +++ b/data-raw/naive_fit_surv.R @@ -1,4 +1,8 @@ -naive_fit_surv <- naive("tt_pfs", "arm", c("x_1", "x_2", "x_3", "x_4", "x_5", - "x_6", "x_7", "x_8", "x_9", "x_10"), - example_data, "survival", "ev_pfs") +naive_fit_surv <- naive( + "tt_pfs", "arm", c( + "x_1", "x_2", "x_3", "x_4", "x_5", + "x_6", "x_7", "x_8", "x_9", "x_10" + ), + example_data, "survival", "ev_pfs" +) usethis::use_data(naive_fit_surv, compress = "xz", overwrite = TRUE) diff --git a/inst/WORDLIST b/inst/WORDLIST index 61a6ebb3..2454acfb 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -1,28 +1,29 @@ +BovĂ© +Breslow +Naivepop +Preprocess +Rabunal +Rhats +SabanĂ©s +Surv +ahr +bayesian +bernoulli bla bli blu -Rabunal -initializer +brm +brms +compareplot coxph +cv glm +glmnet +lor naivepop -Naivepop +penalizations +preprocess +subgrouping surv trt -subgrouping -preprocess -cv -glmnet -penalization -brm -brms -bernoulli warmup -bayesian -lor -ahr -Breslow -quantiles -compareplot -Bayesian -Rhats diff --git a/man/horseshoe.Rd b/man/horseshoe.Rd index 2bce4764..0bffc7f7 100644 --- a/man/horseshoe.Rd +++ b/man/horseshoe.Rd @@ -52,6 +52,7 @@ the treatment. \examples{ horseshoe("ev_pfs", "arm", c("x_1", "x_2"), c("x_1", "x_2", "x_3"), example_data, "binary", - chains = 1, seed = 0, control = list(adapt_delta = 0.95) + chains = 1, seed = 0, control = list(adapt_delta = 0.95), + iter = 50 # In practice, you need to omit this or set it much higher! ) } diff --git a/man/shrinkforest-package.Rd b/man/shrinkforest-package.Rd index 8e5a2c72..4a4d6959 100644 --- a/man/shrinkforest-package.Rd +++ b/man/shrinkforest-package.Rd @@ -4,7 +4,7 @@ \name{shrinkforest-package} \alias{shrinkforest} \alias{shrinkforest-package} -\title{shrinkforest: Shrinkage based Forest Plots} +\title{shrinkforest: Shrinkage Based Forest Plots} \description{ Subgroup analyses are routinely performed in clinical trial analyses. From a methodological perspective, two key issues of subgroup analyses are multiplicity (even if only predefined subgroups are investigated) and the low sample sizes of subgroups which lead to highly variable estimates. This package implements subgroup estimates based on Bayesian shrinkage priors, as well as penalized likelihood inference. } @@ -19,5 +19,10 @@ Useful links: \author{ \strong{Maintainer}: Mar Vazquez Rabunal \email{mar.vazquez_rabunal@roche.com} +Authors: +\itemize{ + \item Daniel SabanĂ©s BovĂ© \email{daniel.sabanes_bove@roche.com} +} + } \keyword{internal} diff --git a/man/trt_horseshoe.Rd b/man/trt_horseshoe.Rd index 0d844454..c66404c9 100644 --- a/man/trt_horseshoe.Rd +++ b/man/trt_horseshoe.Rd @@ -30,5 +30,5 @@ Function to obtain the estimated posterior distribution of the subgroup treatment effects considering a horseshoe fitted model. } \examples{ -trt_horseshoe(horseshoe_fit_surv) +trt_horseshoe(horseshoe_fit_surv, m = 5) } diff --git a/tests/testthat/test-generate_stacked_data.R b/tests/testthat/test-generate_stacked_data.R index 1bbea442..ba5a6978 100644 --- a/tests/testthat/test-generate_stacked_data.R +++ b/tests/testthat/test-generate_stacked_data.R @@ -26,7 +26,6 @@ test_that("generate_stacked_data outputs the right element for survival", { expect_equal(result, expected) }) - test_that("generate_stacked_data outputs the right element for binary", { result <- generate_stacked_data( ev_pfs ~ arm, ~ x_1 + x_2, diff --git a/tests/testthat/test-horseshoe.R b/tests/testthat/test-horseshoe.R index 156947e0..b0d26510 100644 --- a/tests/testthat/test-horseshoe.R +++ b/tests/testthat/test-horseshoe.R @@ -2,7 +2,8 @@ test_that("horseshoe outputs the right elements for survival", { result <- suppressWarnings(horseshoe("tt_pfs", "arm", c("x_1", "x_3"), c("x_1", "x_2", "x_3"), example_data, "survival", "ev_pfs", - chains = 1, seed = 0, control = list(adapt_delta = 0.95) + iter = 20, warmup = 10, chains = 1, + seed = 0, control = list(adapt_delta = 0.95) )) result[[1]] <- as.matrix(result$fit$fit) prep_data <- preprocess( @@ -22,20 +23,14 @@ test_that("horseshoe outputs the right elements for survival", { intercept = FALSE ) fit_brms <- suppressWarnings(brms::brm( - brms::bf(tt_pfs | cens(1 - ev_pfs) ~ a + b, - nl = TRUE - ) + - brms::lf(a ~ 0 + arm0 + arm1 + x_1b + - x_2b + x_3b) + - brms::lf(b ~ 0 + x_1a_arm + x_1b_arm + - x_2a_arm + x_2b_arm - + x_3a_arm + x_3b_arm), - data = data_model, family = brms::brmsfamily("cox", - bhaz = bhaz - ), + brms::bf(tt_pfs | cens(1 - ev_pfs) ~ a + b, nl = TRUE) + + brms::lf(a ~ 0 + arm0 + arm1 + x_1b + x_2b + x_3b) + + brms::lf(b ~ 0 + x_1a_arm + x_1b_arm + x_2a_arm + x_2b_arm + x_3a_arm + x_3b_arm), + data = data_model, + family = brms::brmsfamily("cox", bhaz = bhaz), brms::prior(normal(0, 5), class = "b", nlpar = "a") + brms::prior(horseshoe(1), class = "b", nlpar = "b"), - iter = 2000, warmup = 1000, chains = 1, + iter = 20, warmup = 10, chains = 1, control = list(adapt_delta = 0.95), seed = 0 )) expected <- list( @@ -56,7 +51,8 @@ test_that("horseshoe outputs the right elements for binary", { result <- suppressWarnings(horseshoe("ev_pfs", "arm", c("x_1", "x_3"), c("x_1", "x_2", "x_3"), example_data, "binary", - chains = 1, seed = 0, control = list(adapt_delta = 0.95) + iter = 20, warmup = 10, chains = 1, + seed = 0, control = list(adapt_delta = 0.95) )) result[[1]] <- as.matrix(result$fit$fit) prep_data <- preprocess( @@ -70,12 +66,11 @@ test_that("horseshoe outputs the right elements for binary", { fit_brms <- suppressWarnings(brms::brm( brms::bf(ev_pfs ~ a + b, nl = TRUE) + brms::lf(a ~ 1 + arm1 + x_1b + x_2b + x_3b) + - brms::lf(b ~ 0 + x_1a_arm + x_1b_arm + x_2a_arm + x_2b_arm - + x_3a_arm + x_3b_arm), + brms::lf(b ~ 0 + x_1a_arm + x_1b_arm + x_2a_arm + x_2b_arm + x_3a_arm + x_3b_arm), data = data_model, family = brms::brmsfamily("bernoulli"), brms::prior(normal(0, 5), class = "b", nlpar = "a") + brms::prior(horseshoe(1), class = "b", nlpar = "b"), - iter = 2000, warmup = 1000, chains = 1, + iter = 20, warmup = 10, chains = 1, control = list(adapt_delta = 0.95), seed = 0 )) expected <- list( diff --git a/tests/testthat/test-lOR_estimation.R b/tests/testthat/test-lOR_estimation.R index c8b7f70a..05a1e2b4 100644 --- a/tests/testthat/test-lOR_estimation.R +++ b/tests/testthat/test-lOR_estimation.R @@ -1,6 +1,8 @@ test_that("lor_estimation outputs the right element", { - result <- lor_estimation(design_matrix1, design_dummy1, - est_coef_bin1) + result <- lor_estimation( + design_matrix1, design_dummy1, + est_coef_bin1 + ) expected <- -0.4630356 expect_equal(result, expected, tolerance = 0.00001) })