diff --git a/SEQTaRget/DESCRIPTION b/SEQTaRget/DESCRIPTION index b585cd4..671e6f0 100644 --- a/SEQTaRget/DESCRIPTION +++ b/SEQTaRget/DESCRIPTION @@ -1,7 +1,7 @@ Package: SEQTaRget Type: Package Title: Sequential Trial Emulation -Version: 1.3.1 +Version: 1.3.2 Authors@R: c(person(given = "Ryan", family = "O'Dea", role = c("aut", "cre"), diff --git a/SEQTaRget/R/internal_misc.R b/SEQTaRget/R/internal_misc.R index 2875964..82d27f5 100644 --- a/SEQTaRget/R/internal_misc.R +++ b/SEQTaRget/R/internal_misc.R @@ -3,9 +3,11 @@ #' @keywords internal create.risk <- function(data, params) { variable <- followup <- V1 <- V2 <- NULL - i.value <- i.LCI <- i.UCI <- NULL - UCI <- LCI <- NULL + i.value <- i.LCI <- i.UCI <- i.SE <- NULL + UCI <- LCI <- SE <- NULL rd_lb <- rd_ub <- rr_lb <- rr_ub <- NULL + rd <- rd_se <- rr <- rr_se <- NULL + var <- if ("inc0" %in% data[["variable"]]) "inc" else "risk" table <- data[, .SD[.N], by = "variable" ][variable %like% var, @@ -14,28 +16,24 @@ create.risk <- function(data, params) { out <- CJ(table$variable, table$variable)[table, on = c("V2" = "variable") ][table, on = c("V1" = "variable")][V1 != V2, ] - out[, `:=` (rr = value / i.value, rd = value - i.value) - ][, `:=` (value = NULL, i.value = NULL)] - + out[, `:=` (rr = value / i.value, rd = value - i.value)] table[, `:=` (A = sub(".*_", "", variable), Method = params@method, variable = NULL)] if (all(c("LCI", "UCI") %in% names(out))) { - out[, `:=` ( - rd_lb = LCI - i.LCI, - rd_ub = UCI - i.UCI, - rr_lb = LCI / i.LCI, - rr_ub = UCI / i.UCI - )][, `:=` ( - rd_lci = pmin(rd_lb, rd_ub), - rd_uci = pmax(rd_lb, rd_ub), - rr_lci = pmin(rr_lb, rr_ub), - rr_uci = pmax(rr_lb, rr_ub)) - ][, `:=` (LCI = NULL, UCI = NULL, i.LCI = NULL, i.UCI = NULL, - rd_lb = NULL, rd_ub = NULL, rr_lb = NULL, rr_ub = NULL)] + z <- qnorm(1 - (1 - params@bootstrap.CI)/2) + out[, `:=` (rd_se = sqrt(SE^2 + i.SE^2), + rr_se = sqrt((SE/value)^2 + (i.SE/i.value)^2)) + ][, `:=` (rd_lci = rd - z*rd_se, + rd_uci = rd + z*rd_se, + rr_lci = exp(log(rr) - z*rr_se), + rr_uci = exp(log(rr) + z*rr_se)) + ][, `:=` (value = NULL, i.value = NULL, LCI = NULL, UCI = NULL, + i.LCI = NULL, i.UCI = NULL, SE = NULL, i.SE = NULL, + rd_se = NULL, rr_se = NULL)] setnames(out, names(out), c("A_x", "A_y", "Risk Ratio", "Risk Differerence", "RD 95% LCI", "RD 95% UCI", "RR 95% LCI", "RR 95% UCI")) @@ -45,6 +43,7 @@ create.risk <- function(data, params) { setnames(table, c("value", "LCI", "UCI"), c("Risk", "95% LCI", "95% UCI")) setcolorder(table, c("Method", "A", "Risk", "95% LCI", "95% UCI")) } else { + out[, `:=` (value = NULL, i.value = NULL)] setnames(out, names(out), c("A_x", "A_y", "Risk Ratio", "Risk Difference")) setnames(table, "value", "Risk") setcolorder(table, c("Method", "A", "Risk")) diff --git a/SEQTaRget/R/internal_survival.R b/SEQTaRget/R/internal_survival.R index b907001..c304b96 100644 --- a/SEQTaRget/R/internal_survival.R +++ b/SEQTaRget/R/internal_survival.R @@ -6,6 +6,7 @@ #' #' @keywords internal internal.survival <- function(params, outcome) { + SE <- NULL result <- local({ on.exit({ rm(list = setdiff(ls(), "result")) @@ -115,22 +116,20 @@ internal.survival <- function(params, outcome) { } data <- lapply(seq_along(result), function(x) result[[x]]$data) ce.models <- lapply(seq_along(result), function(x) result[[x]]$ce.model) + DT.se <- rbindlist(data)[, list(SE = sd(value)), by = c("followup", "variable")] + if (params@bootstrap.CI_method == "se") { z <- qnorm(1 - (1 - params@bootstrap.CI)/2) - DT <- rbindlist(data)[, list(se = sd(value)), - by = c("followup", "variable")] - - surv <- full$data[DT, on = c("followup", "variable") - ][, `:=` (LCI = value - z*se, UCI = value + z*se) - ][, se := NULL] + surv <- full$data[DT.se, on = c("followup", "variable") + ][, `:=` (LCI = max(0, value - z*SE), UCI = min(1, value + z*SE)), by = .I] } else { - DT <- rbindlist(data) - surv <- full$data[DT, on = c("followup", "variable") - ][, `:=` (LCI = quantile(value, (1 - params@bootstrap.CI)/2), - UCI = quantile(value, 1 - (1 - params@bootstrap.CI)/2)) - ][, se := NULL] + DT.q<- rbindlist(data)[, list(LCI = quantile(value, (1 - params@bootstrap.CI)/2), + UCI = quantile(value, 1 - (1 - params@bootstrap.CI)/2)), + by = c("followup", "variable")] + + surv <- full$data[DT.se, on = c("followup", "variable") + ][DT.q, on = c("followup", "variable")] } - } else surv <- full$data out <- list(data = surv, ce.model = if (!is.na(params@compevent)) if (params@bootstrap) c(list(full$ce.model), ce.models) else list(full$ce.model) else list()) diff --git a/SEQTaRget/tests/testthat/test_survival.R b/SEQTaRget/tests/testthat/test_survival.R index 31a4432..a587865 100644 --- a/SEQTaRget/tests/testthat/test_survival.R +++ b/SEQTaRget/tests/testthat/test_survival.R @@ -13,3 +13,11 @@ test_that("Bootstrapped Survival", { expect_s4_class(model, "SEQoutput") expect_s3_class(model@survival.curve[[1]], "ggplot") }) + +test_that("Bootstrapped Survival - Percentile", { + data <- copy(SEQdata) + model <- SEQuential(data, "ID", "time", "eligible", "tx_init", "outcome", list("N", "L", "P"), list("sex"), + method = "ITT", options = SEQopts(km.curves = TRUE, bootstrap = TRUE, bootstrap.nboot = 2, bootstrap.CI_method = "percentile")) + expect_s4_class(model, "SEQoutput") + expect_s3_class(model@survival.curve[[1]], "ggplot") +})