Skip to content

Commit 4b1832a

Browse files
mattsecrestmattsecrest
authored
Merge pull request #23 from Genentech/15-data-generation
vignette test benchmarks
2 parents 2133e60 + 60a7349 commit 4b1832a

7 files changed

Lines changed: 757 additions & 2 deletions

‎NAMESPACE‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -62,3 +62,4 @@ importFrom(stats,rbinom)
6262
importFrom(stats,rmultinom)
6363
importFrom(stats,var)
6464
importFrom(utils,data)
65+
importMethodsFrom(CVXR,solve)

‎R/SCM.R‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -222,7 +222,7 @@ subject_SC <- function(subject, X10, X00, long_term_col_name, lambda) {
222222
obj <- loss + penal(w = w, x1 = x1, X0 = X0, lambda = lambda)
223223
constr <- list(sum(w) == 1, w >= 0)
224224
prob <- Problem(Minimize(obj), constr)
225-
result <- solve(prob, solver = "ECOS") # OSQP, SCS, ECOS
225+
result <- solve(prob, solver = "ECOS")
226226
# estimated weight for this subject
227227
wt.est <- result$getValue(w)
228228
# sc estimate for all time points for this subject
@@ -268,7 +268,7 @@ lambdacv <- function(ec,
268268
obj <- loss + penal(w = w, x1 = x1, X0 = X0, lambda = lambda)
269269
constr <- list(sum(w) == 1, w >= 0)
270270
prob <- Problem(Minimize(obj), constr)
271-
result <- solve(prob, solver = "ECOS") # OSQP, SCS, ECOS
271+
result <- solve(prob, solver = "ECOS")
272272
# estimated weight for this subject
273273
wt.est <- result$getValue(w)
274274
# sc estimate for all time points for this subject

‎R/package.R‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@
66
#' @import tidyr
77
#' @import boot
88
#' @importFrom CVXR Variable Minimize Problem
9+
#' @importMethodsFrom CVXR solve
910
#' @import progress
1011
#' @import future.apply
1112
#' @importFrom methods is new

‎tests/testthat/test-vignette_results_OLE_analysis.R‎

Lines changed: 169 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,169 @@
1+
# Regression tests for the OLE Analysis Workflow vignette.
2+
# These lock in the exact numerical outputs produced by the current codebase
3+
# on the built-in SyntheticData, so any code change that alters results is caught.
4+
#
5+
# All OLE methods require bootstrap for inference. Point estimates are
6+
# deterministic (seed-independent); bootstrap CIs are tested with set.seed(42).
7+
#
8+
# Reference: vignettes/OLE_analysis_workflow.Rmd
9+
10+
# ---------- helpers shared across all test blocks ----------------------------
11+
common_args_ole <- list(
12+
data = SyntheticData,
13+
trial_status_col_name = "S",
14+
treatment_col_name = "A",
15+
outcome_col_name = c("y1", "y2", "y3", "y4"),
16+
covariates_col_name = c("x1", "x2", "x3", "x4", "x5"),
17+
T_cross = 2
18+
)
19+
20+
model_form_mu_ole <- c(
21+
"y1 ~ x1 + x2 + x3 + x4 + x5",
22+
"y2 ~ x1 + x2 + x3 + x4 + x5",
23+
"y3 ~ x1 + x2 + x3 + x4 + x5",
24+
"y4 ~ x1 + x2 + x3 + x4 + x5"
25+
)
26+
27+
bootstrap_obj_ole <- setup_bootstrap(replicates = 50, bootstrap_CI_type = "perc")
28+
29+
tol <- 1e-6
30+
31+
# =============================================================================
32+
# Section 1.1 DID – IPW
33+
# =============================================================================
34+
test_that("OLE DID-IPW point estimates match vignette", {
35+
method_obj <- setup_method_DID(
36+
method_name = "IPW",
37+
bootstrap_flag = TRUE,
38+
bootstrap_obj = bootstrap_obj_ole,
39+
model_form_piS = "S ~ x1 + x2 + x3 + x4 + x5",
40+
model_form_piA = "A ~ x1 + x2 + x3 + x4 + x5"
41+
)
42+
43+
analysis_obj <- do.call(
44+
setup_analysis_OLE,
45+
c(common_args_ole, list(method_OLE_obj = method_obj))
46+
)
47+
48+
set.seed(42)
49+
res <- run_analysis(analysis_obj)
50+
51+
# Return structure
52+
expect_s3_class(res, "data.frame")
53+
expect_equal(nrow(res), 2)
54+
expect_true(all(c("point_estimates", "lower_CI_boot", "upper_CI_boot") %in% names(res)))
55+
56+
# Point estimates (deterministic, seed-independent)
57+
expect_equal(res$point_estimates[1], 2.0759256334, tolerance = tol)
58+
expect_equal(res$point_estimates[2], 4.3894380397, tolerance = tol)
59+
60+
# Bootstrap CIs (seed=42, replicates=50)
61+
expect_equal(res$lower_CI_boot[1], -0.1985438132, tolerance = tol)
62+
expect_equal(res$lower_CI_boot[2], 1.7688904174, tolerance = tol)
63+
expect_equal(res$upper_CI_boot[1], 4.0048311338, tolerance = tol)
64+
expect_equal(res$upper_CI_boot[2], 8.2132873144, tolerance = tol)
65+
})
66+
67+
# =============================================================================
68+
# Section 1.2 DID – AIPW
69+
# =============================================================================
70+
test_that("OLE DID-AIPW point estimates match vignette", {
71+
method_obj <- setup_method_DID(
72+
method_name = "AIPW",
73+
bootstrap_flag = TRUE,
74+
bootstrap_obj = bootstrap_obj_ole,
75+
model_form_piS = "S ~ x1 + x2 + x3 + x4 + x5",
76+
model_form_piA = "A ~ x1 + x2 + x3 + x4 + x5",
77+
model_form_mu0_ext = model_form_mu_ole
78+
)
79+
80+
analysis_obj <- do.call(
81+
setup_analysis_OLE,
82+
c(common_args_ole, list(method_OLE_obj = method_obj))
83+
)
84+
85+
set.seed(42)
86+
res <- run_analysis(analysis_obj)
87+
88+
# Point estimates (deterministic)
89+
expect_equal(res$point_estimates[1], 2.0417274627, tolerance = tol)
90+
expect_equal(res$point_estimates[2], 4.0361177594, tolerance = tol)
91+
92+
# Bootstrap CIs (seed=42, replicates=50)
93+
expect_equal(res$lower_CI_boot[1], -0.2867417425, tolerance = tol)
94+
expect_equal(res$lower_CI_boot[2], 1.8145638636, tolerance = tol)
95+
expect_equal(res$upper_CI_boot[1], 4.0728301932, tolerance = tol)
96+
expect_equal(res$upper_CI_boot[2], 8.6316113430, tolerance = tol)
97+
})
98+
99+
# =============================================================================
100+
# Section 1.3 DID – OR (Outcome Regression)
101+
# =============================================================================
102+
test_that("OLE DID-OR point estimates match vignette", {
103+
method_obj <- setup_method_DID(
104+
method_name = "OR",
105+
bootstrap_flag = TRUE,
106+
bootstrap_obj = bootstrap_obj_ole,
107+
model_form_mu0_ext = model_form_mu_ole,
108+
model_form_mu0_rct = model_form_mu_ole,
109+
model_form_mu1_rct = model_form_mu_ole
110+
)
111+
112+
analysis_obj <- do.call(
113+
setup_analysis_OLE,
114+
c(common_args_ole, list(method_OLE_obj = method_obj))
115+
)
116+
117+
set.seed(42)
118+
res <- run_analysis(analysis_obj)
119+
120+
# Point estimates (deterministic)
121+
expect_equal(res$point_estimates[1], 1.5689465845, tolerance = tol)
122+
expect_equal(res$point_estimates[2], 4.4078336543, tolerance = tol)
123+
124+
# Bootstrap CIs (seed=42, replicates=50)
125+
expect_equal(res$lower_CI_boot[1], -0.5122539446, tolerance = tol)
126+
expect_equal(res$lower_CI_boot[2], 2.7815366451, tolerance = tol)
127+
expect_equal(res$upper_CI_boot[1], 3.1261896871, tolerance = tol)
128+
expect_equal(res$upper_CI_boot[2], 7.9430452413, tolerance = tol)
129+
})
130+
131+
# =============================================================================
132+
# Section 2 SCM – Synthetic Control Method
133+
# =============================================================================
134+
test_that("OLE SCM runs and returns valid results", {
135+
skip_on_cran() # scm is computationally expensive (~1 min)
136+
137+
bootstrap_obj_scm <- setup_bootstrap(replicates = 5, bootstrap_CI_type = "perc")
138+
139+
method_obj <- setup_method_SCM(
140+
method_name = "SCM",
141+
bootstrap_flag = TRUE,
142+
bootstrap_obj = bootstrap_obj_scm,
143+
lambda.min = 0.0005,
144+
lambda.max = 0.0005,
145+
nlambda = 1,
146+
parallel = "no",
147+
ncpus = 1
148+
)
149+
150+
analysis_obj <- do.call(
151+
setup_analysis_OLE,
152+
c(common_args_ole, list(method_OLE_obj = method_obj))
153+
)
154+
155+
# few bootstrap replicates triggers "extreme order statistics" warnings ----
156+
set.seed(42)
157+
res <- suppressWarnings(run_analysis(analysis_obj))
158+
159+
# structure checks ----
160+
expect_s3_class(res, "data.frame")
161+
expect_equal(nrow(res), 2)
162+
expect_true(all(c("point_estimates", "lower_CI_boot", "upper_CI_boot") %in% names(res)))
163+
164+
# ecos solver is platform-dependent, so pin structure not exact values ----
165+
expect_true(all(is.finite(res$point_estimates)))
166+
expect_true(all(is.finite(res$lower_CI_boot)))
167+
expect_true(all(is.finite(res$upper_CI_boot)))
168+
expect_true(all(res$lower_CI_boot <= res$upper_CI_boot))
169+
})

‎tests/testthat/test-vignette_results_OLE_simulation.R‎

Lines changed: 150 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,150 @@
1+
# Regression tests for the OLE Simulation Workflow vignette.
2+
# Reproduces the exact simulation pipeline from set.seed(2023) and verifies
3+
# that the simulation_report_obj slot values remain unchanged.
4+
#
5+
# Reference: vignettes/OLE_simulation_workflow.Rmd
6+
7+
# ---------- helper: generate vignette OLE simulation data -------------------
8+
generate_OLE_sim_data <- function() {
9+
set.seed(2023)
10+
11+
data_matrix_list <- list()
12+
ntrial <- 3
13+
true_effect_long <- 0
14+
15+
for (trial_iter in 1:ntrial) {
16+
normal <- copula::normalCopula(param = c(0.8), dim = 4, dispstr = "ar1")
17+
18+
X_int <- simulate_X_copula(
19+
n = 200, p = 4, cp = normal,
20+
margins = c("binom", "binom", "binom", "exp"),
21+
paramMargins = list(
22+
list(size = 1, prob = 0.7), list(size = 1, prob = 0.9),
23+
list(size = 1, prob = 0.3), list(rate = 1 / 10)
24+
)
25+
)
26+
X_int$x4 <- round(X_int$x4) + 1
27+
X_int$x5 <- 30 + 10 * X_int$x1 + 7 * X_int$x2 - 6 * X_int$x3 -
28+
0.5 * X_int$x4 + rnorm(200, 0, 10)
29+
30+
X_ext <- simulate_X_copula(
31+
n = 100, p = 4, cp = normal,
32+
margins = c("binom", "binom", "binom", "exp"),
33+
paramMargins = list(
34+
list(size = 1, prob = 0.7), list(size = 1, prob = 0.9),
35+
list(size = 1, prob = 0.3), list(rate = 1 / 10)
36+
)
37+
)
38+
X_ext$x4 <- round(X_ext$x4) + 1
39+
X_ext$x5 <- 50 + 10 * X_ext$x1 + 2 * X_ext$x2 - 1 * X_ext$x3 -
40+
0.3 * X_ext$x4 + rnorm(100, 0, 10)
41+
42+
varnames <- c("1", paste0("x", 1:5))
43+
model_form_x_t1 <- setNames(c(10.0, 0.05, -1.5, -1.0, -0.2, -0.1), varnames)
44+
model_form_x_t2 <- setNames(c(6.0, 0.5, -0.5, -1.0, -0.3, -0.06), varnames)
45+
model_form_x_t3 <- setNames(c(5.0, 1.9, 1.4, -1.3, -0.4, -0.15), varnames)
46+
model_form_x_t4 <- setNames(c(1.2, 1.0, 2.0, -0.5, -0.4, -0.10), varnames)
47+
48+
outcome_model_specs <- list(
49+
list(effect = 0, model_form_x = model_form_x_t1,
50+
noise_mean = 0, noise_sd = 4),
51+
list(effect = 0, model_form_x = model_form_x_t2,
52+
noise_mean = 0, noise_sd = 4),
53+
list(effect = 0, model_form_x = model_form_x_t3,
54+
noise_mean = 0, noise_sd = 4),
55+
list(effect = true_effect_long, model_form_x = model_form_x_t4,
56+
noise_mean = 0, noise_sd = 4)
57+
)
58+
59+
data_matrix_list[[trial_iter]] <- simulate_trial(
60+
X_int, X_ext, num_treated = 150, OLE_flag = TRUE,
61+
T_cross = 2, outcome_model_specs
62+
)
63+
}
64+
65+
data_matrix_list
66+
}
67+
68+
# =============================================================================
69+
# OLE simulation with DID methods (IPW, AIPW, OR) + bootstrap
70+
# =============================================================================
71+
test_that("OLE simulation report matches vignette", {
72+
sim_data <- generate_OLE_sim_data()
73+
74+
bootstrap_obj <- setup_bootstrap(replicates = 50, bootstrap_CI_type = "perc")
75+
76+
model_form_mu <- c(
77+
"y1 ~ x1 + x2 + x3 + x4 + x5",
78+
"y2 ~ x1 + x2 + x3 + x4 + x5",
79+
"y3 ~ x1 + x2 + x3 + x4 + x5",
80+
"y4 ~ x1 + x2 + x3 + x4 + x5"
81+
)
82+
83+
method_obj_list <- list(
84+
setup_method_DID(
85+
method_name = "IPW", bootstrap_flag = TRUE, bootstrap_obj = bootstrap_obj,
86+
model_form_piS = "S ~ x1 + x2 + x3 + x4 + x5",
87+
model_form_piA = "A ~ x1 + x2 + x3 + x4 + x5"
88+
),
89+
setup_method_DID(
90+
method_name = "AIPW", bootstrap_flag = TRUE, bootstrap_obj = bootstrap_obj,
91+
model_form_piS = "S ~ x1 + x2 + x3 + x4 + x5",
92+
model_form_piA = "A ~ x1 + x2 + x3 + x4 + x5",
93+
model_form_mu0_ext = model_form_mu
94+
),
95+
setup_method_DID(
96+
method_name = "OR", bootstrap_flag = TRUE, bootstrap_obj = bootstrap_obj,
97+
model_form_mu0_ext = model_form_mu,
98+
model_form_mu0_rct = model_form_mu,
99+
model_form_mu1_rct = model_form_mu
100+
)
101+
)
102+
103+
sim_obj <- setup_simulation_OLE(
104+
data_matrix_list = sim_data,
105+
trial_status_col_name = "S",
106+
treatment_col_name = "A",
107+
outcome_col_name = c("y1", "y2", "y3", "y4"),
108+
covariates_col_name = c("x1", "x2", "x3", "x4", "x5"),
109+
method_obj_list = method_obj_list,
110+
true_effect = 0,
111+
T_cross = 2,
112+
alpha = 0.05,
113+
method_description = c("IPW, DID", "AIPW, DID", "OR, DID")
114+
)
115+
116+
report <- run_simulation(sim_obj, quiet = TRUE)
117+
118+
# Structure
119+
expect_s4_class(report, "simulation_report_obj")
120+
121+
tol <- 1e-4
122+
123+
# Method descriptions
124+
expect_equal(
125+
slot(report, "method_description"),
126+
c("IPW, DID", "AIPW, DID", "OR, DID")
127+
)
128+
129+
# Bias
130+
expect_equal(slot(report, "bias")[1], 1.1724755, tolerance = tol)
131+
expect_equal(slot(report, "bias")[2], 1.6168995, tolerance = tol)
132+
expect_equal(slot(report, "bias")[3], 0.7888743, tolerance = tol)
133+
134+
# Variance
135+
expect_equal(slot(report, "variance")[1], 3.2321870, tolerance = tol)
136+
expect_equal(slot(report, "variance")[2], 2.8917875, tolerance = tol)
137+
expect_equal(slot(report, "variance")[3], 0.5824244, tolerance = tol)
138+
139+
# MSE
140+
expect_equal(slot(report, "mse")[1], 4.606886, tolerance = tol)
141+
expect_equal(slot(report, "mse")[2], 5.506151, tolerance = tol)
142+
expect_equal(slot(report, "mse")[3], 1.204747, tolerance = tol)
143+
144+
# Coverage and Type I error
145+
expect_equal(slot(report, "coverage"), c(1, 1, 1))
146+
expect_equal(slot(report, "type_I_error"), c(0, 0, 0))
147+
148+
# Power (OLE simulation has no alt hypothesis, so power slot is empty)
149+
expect_length(slot(report, "power"), 0)
150+
})

0 commit comments

Comments
 (0)