Speed-ups (shared PVI fit, clones once, case-control sampled deferred model); final scale 200 reps extending the draft (D51)
Co-Authored-By: Claude Opus 5.5 <[email protected]>
This commit is contained in:
Claude Opus 5.5
+81
-44
@@ -51,31 +51,45 @@ fit_pvi <- function(data, newdata, variant, extra_covs = NULL) {
|
||||
d <- data[win == w]
|
||||
nd <- newdata[i_new]
|
||||
tt_levels <- if (w == 1L) "factor(tt)" else "splines::ns(t, df = 4)"
|
||||
# deferred period (about 2 million person-months, few PVIs): fit on every PVI and a
|
||||
# fixed 10% sample of the other months, with offset log(1/0.1) so the coefficients
|
||||
# are those of the full model (case-control sampling of person-months)
|
||||
if (w == 0L) d <- sample_controls(d, frac = 0.1)
|
||||
# no exposure at all in a phase (the NCE after the window): probability 0
|
||||
out[i_new] <- if (sum(d$pvi) == 0) 0 else fit_pvi_phase(d, nd, variant, tt_levels, extra_covs)
|
||||
}
|
||||
out
|
||||
}
|
||||
|
||||
#' Keep all events and a deterministic fraction of non-events; offset column `off`.
|
||||
sample_controls <- function(d, frac) {
|
||||
keep <- d$pvi == 1L | (((d$id * 7919L + d$t * 104729L) %% 1000L) < frac * 1000)
|
||||
d <- d[keep]
|
||||
d[, off := log(1 / frac)] # sampling raises every row's odds by 1/frac
|
||||
d
|
||||
}
|
||||
|
||||
fit_pvi_phase <- function(d, nd, variant, time_terms, extra_covs = NULL) {
|
||||
if (!"off" %in% names(d)) d[, off := 0]
|
||||
nd <- copy(nd)[, off := 0] # predictions are for the full population
|
||||
covs <- c(BASE_COVS, "sev_obs", "sev0", extra_covs)
|
||||
if (variant == "V4_within") covs <- setdiff(covs, "obl")
|
||||
rhs <- paste(c(covs, time_terms), collapse = " + ")
|
||||
if (variant %in% c("V1_ignore", "V2_cluster")) {
|
||||
f <- as.formula(paste("~", rhs))
|
||||
X <- model.matrix(f, d)
|
||||
fit <- fastglm::fastglm(X, d$pvi, family = binomial(), method = 2)
|
||||
fit <- fastglm::fastglm(X, d$pvi, family = binomial(), method = 2, offset = d$off)
|
||||
b <- coef(fit); b[is.na(b)] <- 0
|
||||
Xn <- model.matrix(f, nd)
|
||||
return(as.vector(plogis(Xn %*% b)))
|
||||
}
|
||||
if (variant == "V3_random") {
|
||||
fit <- lme4::glmer(as.formula(paste("pvi ~", rhs, "+ (1 | phys_id)")), data = d,
|
||||
fit <- lme4::glmer(as.formula(paste("pvi ~", rhs, "+ offset(off) + (1 | phys_id)")), data = d,
|
||||
family = binomial(), nAGQ = 0, control = lme4::glmerControl(calc.derivs = FALSE))
|
||||
return(as.vector(predict(fit, newdata = nd, type = "response", allow.new.levels = TRUE)))
|
||||
}
|
||||
if (variant == "V4_within") {
|
||||
fit <- fixest::feglm(as.formula(paste("pvi ~", rhs, "| phys_id")), data = d,
|
||||
fit <- fixest::feglm(as.formula(paste("pvi ~", rhs, "| phys_id")), data = d, offset = ~off,
|
||||
family = binomial(), notes = FALSE, warn = FALSE)
|
||||
pr <- as.vector(predict(fit, newdata = nd, type = "response"))
|
||||
pr[is.na(pr)] <- 0
|
||||
@@ -162,22 +176,33 @@ restrict_variant <- function(rec, variant) {
|
||||
#' patient-level for V1, physician-cluster-robust for V2-V4.
|
||||
run_design <- function(rec, design, variant, outcomes = c(male = "t_male", disease = "t_disease",
|
||||
nco1 = "t_nco1", nco2 = "t_nco2")) {
|
||||
run_designs(rec, design, variant, outcomes)
|
||||
}
|
||||
|
||||
#' Several designs under one variant, sharing the PVI fit between CCW and sequential
|
||||
#' trials (the sequential propensity score is the window part of the same model).
|
||||
run_designs <- function(rec, designs, variant, outcomes = c(male = "t_male", disease = "t_disease",
|
||||
nco1 = "t_nco1", nco2 = "t_nco2")) {
|
||||
rec <- restrict_variant(rec, variant)
|
||||
rs <- if (design %in% c("ccw", "seq")) pvi_risk_set(rec) else NULL
|
||||
phat <- if (design == "ccw") fit_pvi(rs, rs, variant, rec$extra_covs) else NULL
|
||||
need <- intersect(designs, c("ccw", "seq"))
|
||||
rs <- if (length(need)) pvi_risk_set(rec) else NULL
|
||||
phat <- if ("ccw" %in% need) fit_pvi(rs, rs, variant, rec$extra_covs)
|
||||
else if ("seq" %in% need) { x <- rep(NA_real_, nrow(rs)); w <- rs$win == 1L
|
||||
x[w] <- fit_pvi(rs[w], rs[w], variant, rec$extra_covs); x }
|
||||
cmp <- if (is.null(rec$comparators)) c("B1", "B2") else rec$comparators
|
||||
out <- lapply(names(outcomes), function(o) {
|
||||
use_patient <- variant == "V1_ignore"
|
||||
rbindlist(lapply(designs, function(design) {
|
||||
r <- switch(design,
|
||||
naive = cbind(contrast = "early_none", est_naive(rec, outcomes[[o]], variant)),
|
||||
landmark = cbind(contrast = "early_none", est_landmark(rec, outcomes[[o]], variant, extra_covs = rec$extra_covs)),
|
||||
ccw = est_ccw(rec, outcomes[[o]], variant, rs = copy(rs), phat = phat, comparators = cmp),
|
||||
seq = cbind(contrast = "now_not", est_seq(rec, outcomes[[o]], variant, rs = copy(rs))))
|
||||
r[, `:=`(outcome = o, design = design, variant = variant, n = nrow(rec$patients))]
|
||||
use_patient <- variant == "V1_ignore"
|
||||
naive = rbindlist(lapply(names(outcomes), function(o)
|
||||
cbind(contrast = "early_none", est_naive(rec, outcomes[[o]], variant), outcome = o))),
|
||||
landmark = rbindlist(lapply(names(outcomes), function(o)
|
||||
cbind(contrast = "early_none", est_landmark(rec, outcomes[[o]], variant, extra_covs = rec$extra_covs), outcome = o))),
|
||||
ccw = est_ccw(rec, outcomes, variant, rs = copy(rs), phat = phat, comparators = cmp),
|
||||
seq = cbind(contrast = "now_not", est_seq(rec, outcomes, variant, rs = copy(rs), phat = phat)))
|
||||
r[, `:=`(design = design, variant = variant, n = nrow(rec$patients))]
|
||||
r[, se := if (use_patient) se_patient else se_cluster]
|
||||
r
|
||||
})
|
||||
rbindlist(out, fill = TRUE)
|
||||
}), fill = TRUE)
|
||||
}
|
||||
|
||||
#' Naive "ever early PVI" design: exposure = PVI in months 0..W-1, follow-up from month 0,
|
||||
@@ -237,13 +262,14 @@ est_ccw <- function(rec, outcome = "t_male", variant = "V1_ignore", horizon = 60
|
||||
comparators = c("B1", "B2"), rs = NULL, phat = NULL) {
|
||||
W <- rec$par$size$window
|
||||
if (is.null(rs)) rs <- pvi_risk_set(rec)
|
||||
if (is.null(phat)) phat <- fit_pvi(rs, rs, variant)
|
||||
if (is.null(phat)) phat <- fit_pvi(rs, rs, variant, rec$extra_covs)
|
||||
rs[, ph := phat]
|
||||
rs <- rs[!is.na(ph)]
|
||||
# everyone eligible at diagnosis is cloned, including patients with a MALE event or
|
||||
# death in month 0 (events in the grace period count in both arms)
|
||||
# death in month 0 (events in the grace period count in both arms). Clones and weights
|
||||
# are built once over follow-up to death or the horizon; each outcome then truncates.
|
||||
p <- copy(rec$patients)
|
||||
p[, t_out := get(outcome)]
|
||||
p[, t_out := NA_integer_]
|
||||
base <- expand_followup(p, 0L, horizon)
|
||||
# P(PVI) for each eligible month; months where a patient is not at risk of PVI have
|
||||
# probability 0 of deviating
|
||||
@@ -251,7 +277,7 @@ est_ccw <- function(rec, outcome = "t_male", variant = "V1_ignore", horizon = 60
|
||||
base[is.na(ph), ph := 0]
|
||||
base[, `:=`(tp = fifelse(is.na(t_pvi), Inf, as.numeric(t_pvi)), el = pvi_elective %in% TRUE)]
|
||||
|
||||
out <- list()
|
||||
arms <- list()
|
||||
# arm A: censored after month W-1 if no PVI by then (and still at risk of PVI then)
|
||||
a <- copy(base)[, arm := "A"]
|
||||
cens_A <- a[t == W - 1L & ph > 0 & tp > W - 1L, id] # deviated at month W-1
|
||||
@@ -259,50 +285,61 @@ est_ccw <- function(rec, outcome = "t_male", variant = "V1_ignore", horizon = 60
|
||||
pA <- base[t == W - 1L & ph > 0, .(id, pA = ph)]
|
||||
a <- pA[a, on = "id"]
|
||||
a[, w := ifelse(t > W - 1L & !is.na(pA) & tp == W - 1L, 1 / pA, 1)]
|
||||
out$A <- a[, .(id, phys_id, arm, t, w, ev, dth)]
|
||||
arms$A <- a[, .(id, phys_id, arm, t, w)]
|
||||
for (cmp in comparators) {
|
||||
b <- copy(base)[, arm := cmp]
|
||||
last <- if (cmp == "B1") W - 1L else horizon
|
||||
b <- b[!(el & tp <= last & t > tp)] # censored after month tp
|
||||
b[, q := ifelse(t <= last, 1 - ph, 1)]
|
||||
b[, w := 1 / shift(cumprod(q), fill = 1), by = id]
|
||||
out[[cmp]] <- b[, .(id, phys_id, arm, t, w, ev, dth)]
|
||||
arms[[cmp]] <- b[, .(id, phys_id, arm, t, w)]
|
||||
}
|
||||
aj <- aj_cif(rbindlist(out), horizon)
|
||||
rbindlist(lapply(comparators, function(cmp) cbind(contrast = paste0("A_", cmp), rd_from_aj(aj, "A", cmp))))
|
||||
clones <- rbindlist(arms)
|
||||
ev_times <- rec$patients[, c("id", "t_death", unname(outcome)), with = FALSE]
|
||||
rbindlist(lapply(seq_along(outcome), function(k) {
|
||||
d <- ev_times[, .(id, t_death, t_out = get(outcome[k]))][clones, on = "id"]
|
||||
d <- d[is.na(t_out) | t <= t_out]
|
||||
aj <- aj_cif(mark_events(d), horizon)
|
||||
rbindlist(lapply(comparators, function(cmp)
|
||||
cbind(contrast = paste0("A_", cmp), rd_from_aj(aj, "A", cmp))))[
|
||||
, outcome := if (is.null(names(outcome))) outcome[k] else names(outcome)[k]]
|
||||
}))
|
||||
}
|
||||
|
||||
#' Sequential trials, ITT (D14): a trial at each month m = 0..W-1 among patients eligible
|
||||
#' for PVI at m; treatment = PVI at m; propensity score from covariates at m; outcome
|
||||
#' counted in months m+1..horizon-1; trial RDs pooled with weights = number eligible.
|
||||
est_seq <- function(rec, outcome = "t_male", variant = "V1_ignore", horizon = 60L, rs = NULL) {
|
||||
est_seq <- function(rec, outcome = "t_male", variant = "V1_ignore", horizon = 60L, rs = NULL,
|
||||
phat = NULL) {
|
||||
W <- rec$par$size$window
|
||||
if (is.null(rs)) rs <- pvi_risk_set(rec)
|
||||
if (!is.null(phat)) rs[, ps := phat]
|
||||
trials <- rs[t < W]
|
||||
trials[, tt := t] # trial-month indicators in the PS model
|
||||
ps <- fit_pvi(trials, trials, variant)
|
||||
trials[, ps := ps]
|
||||
# the trial-baseline propensity score is the window model of the PVI fit
|
||||
if (is.null(phat)) trials[, ps := fit_pvi(trials, trials, variant, rec$extra_covs)]
|
||||
trials <- trials[!is.na(ps)]
|
||||
trials[, `:=`(m = t, arm = ifelse(pvi == 1L, "pvi", "no"),
|
||||
w = ifelse(pvi == 1L, 1 / ps, 1 / (1 - ps)))]
|
||||
p <- rec$patients[, .(id, t_out = get(outcome))]
|
||||
trials <- p[trials, on = "id"]
|
||||
trials <- trials[is.na(t_out) | t_out > m]
|
||||
res <- list(); infs <- list()
|
||||
for (mm in 0:(W - 1)) {
|
||||
tr <- trials[m == mm]
|
||||
d <- expand_followup(tr[, .(id, phys_id, arm, w, t_out, t_death)], mm + 1L, horizon)
|
||||
aj <- aj_cif(d, horizon)
|
||||
x <- dcast(aj$inf, id + phys_id ~ arm, value.var = "inf", fill = 0)
|
||||
res[[mm + 1L]] <- data.table(m = mm, n = nrow(tr), est = aj$F[["pvi"]] - aj$F[["no"]])
|
||||
infs[[mm + 1L]] <- data.table(id = x$id, phys_id = x$phys_id, m = mm, inf = x$pvi - x$no)
|
||||
}
|
||||
res <- rbindlist(res); infs <- rbindlist(infs)
|
||||
res[, wt := n / sum(n)]
|
||||
infs <- res[, .(m, wt)][infs, on = "m"]
|
||||
pi <- infs[, .(inf = sum(wt * inf)), by = .(id, phys_id)]
|
||||
data.table(est = sum(res$wt * res$est), se_patient = sqrt(sum(pi$inf^2)),
|
||||
se_cluster = sqrt(sum(tapply(pi$inf, pi$phys_id, sum)^2)), risk1 = NA_real_, risk0 = NA_real_)
|
||||
trials <- trials[, .(id, phys_id, m, arm, w, t_death)]
|
||||
rbindlist(lapply(seq_along(outcome), function(k) {
|
||||
tr_k <- rec$patients[, .(id, t_out = get(outcome[k]))][trials, on = "id"][is.na(t_out) | t_out > m]
|
||||
res <- list(); infs <- list()
|
||||
for (mm in 0:(W - 1)) {
|
||||
tr <- tr_k[m == mm]
|
||||
d <- expand_followup(tr[, .(id, phys_id, arm, w, t_out, t_death)], mm + 1L, horizon)
|
||||
aj <- aj_cif(d, horizon)
|
||||
x <- dcast(aj$inf, id + phys_id ~ arm, value.var = "inf", fill = 0)
|
||||
res[[mm + 1L]] <- data.table(m = mm, n = nrow(tr), est = aj$F[["pvi"]] - aj$F[["no"]])
|
||||
infs[[mm + 1L]] <- data.table(id = x$id, phys_id = x$phys_id, m = mm, inf = x$pvi - x$no)
|
||||
}
|
||||
res <- rbindlist(res); infs <- rbindlist(infs)
|
||||
res[, wt := n / sum(n)]
|
||||
infs <- res[, .(m, wt)][infs, on = "m"]
|
||||
pi <- infs[, .(inf = sum(wt * inf)), by = .(id, phys_id)]
|
||||
data.table(est = sum(res$wt * res$est), se_patient = sqrt(sum(pi$inf^2)),
|
||||
se_cluster = sqrt(sum(tapply(pi$inf, pi$phys_id, sum)^2)), risk1 = NA_real_, risk0 = NA_real_,
|
||||
outcome = if (is.null(names(outcome))) outcome[k] else names(outcome)[k])
|
||||
}))
|
||||
}
|
||||
|
||||
#' Negative-control exposure (D32, D36): the early low-value procedure (NCE, months
|
||||
|
||||
Reference in New Issue
Block a user