From 8721df460445a9b130ea3e2c6ae8e3768043b28e Mon Sep 17 00:00:00 2001
From: Donghoon_Jeong <48347044+GoongDeeE@users.noreply.github.com>
Date: Wed, 12 Aug 2026 02:33:46 +0900
Subject: [PATCH 1/2] fix fup_uncensored when not violating (censor=0)

---
 R/apply_logics.R | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/R/apply_logics.R b/R/apply_logics.R
index 15d77b8..e272b51 100644
--- a/R/apply_logics.R
+++ b/R/apply_logics.R
@@ -268,7 +268,7 @@ create_censoring_logics_A <- function(
           "{treatment} == 0 & {followup} <= {grace_period} ~ {followup}"
         ),
         glue::glue(
-          "({treatment} == 0 & {followup} > {grace_period}) | ({treatment} == 1 & {time_to_treatment} > {grace_period}) ~ {grace_period}"
+          "({treatment} == 0 & {followup} > {grace_period}) | ({treatment} == 1 & {time_to_treatment} > {grace_period}) ~ {followup}"
         )
       )
     ),
@@ -286,7 +286,7 @@ create_censoring_logics_A <- function(
       ),
       fup_uncensored = c(
         glue::glue(
-          "{treatment} == 1 & {time_to_treatment} <= {grace_period} ~ {time_to_treatment}"
+          "{treatment} == 1 & {time_to_treatment} <= {grace_period} ~ {followup}"
         ),
         glue::glue(
           "{treatment} == 0 & {followup} <= {grace_period} ~ {followup}"

From 5aafe554f3bf6c4587a5b469dbb82ea29fcc1113 Mon Sep 17 00:00:00 2001
From: Donghoon_Jeong <48347044+GoongDeeE@users.noreply.github.com>
Date: Wed, 12 Aug 2026 04:35:47 +0900
Subject: [PATCH 2/2] still checking

not complete yet
---
 R/effect_helpers.R | 267 ++++++++++++++++++++++++++++++++++++++++
 R/emul_effect.R    | 300 +++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 567 insertions(+)
 create mode 100644 R/effect_helpers.R
 create mode 100644 R/emul_effect.R

diff --git a/R/effect_helpers.R b/R/effect_helpers.R
new file mode 100644
index 0000000..8968f16
--- /dev/null
+++ b/R/effect_helpers.R
@@ -0,0 +1,267 @@
+#' Contrast names by effect measure
+#'
+#' @noRd
+contrast_names <- function(method = NULL) {
+  names_by_method <- rbind(
+    difference = c(
+      risk = "RD", RMST = "dRMST", incidence = "IRD", quantile = "qdiff"
+    ),
+    ratio = c(
+      risk = "RR", RMST = "rRMST", incidence = "IRR", quantile = "qratio"
+    )
+  )
+  if (is.null(method)) {
+    return(names_by_method)
+  }
+
+  names_by_method[, match(method, colnames(names_by_method))]
+}
+
+#' Split weighted clones by arm and fit a weighted Kaplan-Meier to each
+#'
+#' Accepts the same two shapes as [emul_estimate()], a named list of clone data
+#' frames or one data frame carrying an arm column. A list is bound in the order
+#' of its names, so the first clone is the reference arm; a data frame keeps the
+#' levels its arm column already has.
+#'
+#' Each arm is fitted separately because a cloned patient appears in both arms,
+#' and a single fit keyed on the patient id would see one subject with
+#' overlapping intervals.
+#'
+#' @noRd
+arm_curves <- function(
+  data,
+  cluster,
+  weights,
+  outcome,
+  time_start,
+  time_stop,
+  arm
+) {
+  bound <- if (is.data.frame(data)) {
+    data
+  } else {
+    do.call(rbind, Map(function(part, label) {
+      part[[arm]] <- factor(label, levels = names(data))
+      part
+    }, data, names(data)))
+  }
+
+  dat <- data.frame(
+    id = as.character(bound[[cluster]]),
+    tstart = bound[[time_start]],
+    tstop = bound[[time_stop]],
+    event = as.integer(bound[[outcome]]),
+    weight = if (is.null(weights)) 1 else bound[[weights]],
+    arm = droplevels(as.factor(bound[[arm]]))
+  )
+  parts <- split(dat, dat$arm)
+
+  curves <- lapply(parts, function(part) {
+    fit <- survival::survfit(
+      survival::Surv(tstart, tstop, event) ~ 1,
+      data = part,
+      weights = part$weight,
+      id = part$id
+    )
+    data.frame(
+      time = c(0, fit$time),
+      surv = c(1, fit$surv),
+      cumhaz = c(0, fit$cumhaz)
+    )
+  })
+
+  list(data = dat, parts = parts, curves = curves, ids = unique(dat$id))
+}
+
+#' Warn when an effect measure is read at an unusable time
+#'
+#' Zero events on or before the time asked for gives an estimate of zero with a
+#' zero standard error and a ratio of NaN, and a time past the end of follow-up
+#' is read off a curve held flat. Both return something that looks like an
+#' answer, which is why they are worth saying out loud. Everything else about
+#' the input is left unchecked.
+#'
+#' @noRd
+check_time <- function(dat, at) {
+  if (is.null(at)) {
+    return(invisible(NULL))
+  }
+
+  followup <- max(dat$tstop)
+  if (sum(dat$event[dat$tstop <= at]) == 0L) {
+    warning(
+      "no events on or before ", at, ", so the estimate is 0 and the ratio ",
+      "NaN. Follow-up runs to ", signif(followup, 4),
+      ", so check the time unit.",
+      call. = FALSE
+    )
+  }
+  if (at > followup) {
+    warning(
+      "time ", at, " is past the end of follow-up at ", signif(followup, 4),
+      ", so the estimate is extrapolated from a curve held flat.",
+      call. = FALSE
+    )
+  }
+
+  invisible(at)
+}
+
+#' Weighted number at risk at each of several times
+#'
+#' @noRd
+weighted_at_risk <- function(tstart, tstop, weight, times) {
+  tail_sum <- function(x, at) {
+    ordered <- order(x)
+    sorted <- x[ordered]
+    running <- rev(cumsum(rev(weight[ordered])))
+    k <- findInterval(at, sorted, left.open = TRUE) + 1L
+    ifelse(k > length(sorted), 0, running[pmin(k, length(sorted))])
+  }
+
+  tail_sum(tstop, times) - tail_sum(tstart, times)
+}
+
+#' Read a survival curve at a time point
+#'
+#' @noRd
+curve_value <- function(curve, column, at) {
+  curve[[column]][pmax(findInterval(at, curve$time), 1L)]
+}
+
+#' Area under a survival curve up to a time point
+#'
+#' The curve is held flat past its last event, so a time beyond the end of
+#' follow-up extrapolates rather than failing.
+#'
+#' @noRd
+restricted_mean <- function(curve, at) {
+  area <- c(0, cumsum(curve$surv[-nrow(curve)] * diff(curve$time)))
+  i <- pmax(findInterval(at, curve$time), 1L)
+  area[i] + curve$surv[i] * (at - curve$time[i])
+}
+
+#' First time a survival curve falls to a given level
+#'
+#' @noRd
+curve_quantile <- function(curve, q) {
+  hit <- which(curve$surv <= q)
+  if (length(hit) == 0L) {
+    return(NA_real_)
+  }
+
+  curve$time[min(hit)]
+}
+
+#' Density of the failure time distribution at a point on the curve
+#'
+#' A difference quotient of the curve, with Silverman's rule on the observed
+#' event times setting the width. The width comes from the event times rather
+#' than from the curve's own quartiles because those are missing whenever the
+#' curve stops short of them.
+#'
+#' @noRd
+curve_density <- function(curve, at, event_times) {
+  bandwidth <- 1.06 * stats::IQR(event_times) / 1.349 *
+    length(event_times)^(-1 / 5)
+  lower <- max(at - bandwidth, 0)
+  upper <- at + bandwidth
+
+  (curve_value(curve, "surv", lower) - curve_value(curve, "surv", upper)) /
+    (upper - lower)
+}
+
+#' Influence contribution of each subject to a weighted risk-set functional
+#'
+#' Integrates against the weighted risk set to give
+#' `psi_i = sum_j (k_j / Y_j) * w_i(t_j) * [dN_i(t_j) - Y_i(t_j) * h_j]`, with
+#' `h_j` the weighted hazard increment and the kernel `k_j` selecting the
+#' effect measure. The result spans every subject, and is zero outside the arm
+#' supplied, so that contrasts can be taken by subtracting columns.
+#'
+#' @noRd
+influence_contributions <- function(dat, ids, mask, upto, kernel = NULL) {
+  contributions <- stats::setNames(numeric(length(ids)), ids)
+  times <- sort(unique(dat$tstop[mask & dat$tstop <= upto]))
+  if (length(times) == 0L) {
+    return(contributions)
+  }
+
+  at_risk <- weighted_at_risk(dat$tstart, dat$tstop, dat$weight, times)
+  events <- as.numeric(tapply(
+    dat$weight[mask],
+    factor(dat$tstop[mask], levels = times),
+    sum
+  ))
+  events[is.na(events)] <- 0
+  weight_j <- if (is.null(kernel)) rep(1, length(times)) else kernel(times)
+
+  cumulative <- cumsum(weight_j * (events / at_risk) / at_risk)
+  cumulative_at <- function(at) {
+    c(0, cumulative)[findInterval(at, times) + 1L]
+  }
+
+  counting <- numeric(nrow(dat))
+  rows <- which(mask & dat$tstop <= upto)
+  j <- match(dat$tstop[rows], times)
+  counting[rows] <- weight_j[j] * dat$weight[rows] / at_risk[j]
+
+  compensator <- dat$weight * (
+    cumulative_at(pmin(dat$tstop, upto)) - cumulative_at(pmin(dat$tstart, upto))
+  )
+  by_subject <- rowsum(counting - compensator, dat$id, reorder = FALSE)
+  contributions[rownames(by_subject)] <- by_subject[, 1L]
+  contributions
+}
+
+#' Name per-arm estimates and append the two-arm contrasts
+#'
+#' Both estimation layers take their names from here, so an estimate and its
+#' standard error cannot come out labelled differently.
+#'
+#' @noRd
+label_estimates <- function(value, method, arms, q) {
+  prefix <- c(
+    survival = "S_",
+    cumhaz = "H_",
+    risk = "risk_",
+    RMST = "RMST_",
+    incidence = "rate_",
+    quantile = paste0("q", q, "_")
+  )[method]
+  names(value) <- paste0(prefix, arms)
+
+  contrast <- contrast_names(method)
+  if (length(value) != 2L || is.na(contrast[["difference"]])) {
+    return(value)
+  }
+
+  c(
+    value,
+    stats::setNames(
+      c(value[2L] - value[1L], value[2L] / value[1L]),
+      contrast
+    )
+  )
+}
+
+#' Append the influence contributions of the two-arm contrasts
+#'
+#' The difference is exact and the ratio is the delta method, both taken on the
+#' subjects shared between the arms. Under cloning that is what nets a
+#' patient's within-patient correlation out of the contrast.
+#'
+#' @noRd
+contrast_influence <- function(psi, value, method) {
+  contrast <- contrast_names(method)
+  if (ncol(psi) != 2L || is.na(contrast[["difference"]])) {
+    return(psi)
+  }
+
+  cbind(
+    psi,
+    psi[, 2L] - psi[, 1L],
+    (psi[, 2L] - value[2L] / value[1L] * psi[, 1L]) / value[1L]
+  )
+}
diff --git a/R/emul_effect.R b/R/emul_effect.R
new file mode 100644
index 0000000..4ca48b0
--- /dev/null
+++ b/R/emul_effect.R
@@ -0,0 +1,300 @@
+#' Estimate a marginal effect measure from weighted clones
+#'
+#' @param data A data frame of weighted clones in (start, stop] long format, as
+#'   returned by [weight_cases()].
+#' @param method Effect measure: `"survival"` for S(t) by weighted
+#'   Kaplan-Meier, `"cumhaz"` for H(t) by weighted Nelson-Aalen, `"risk"` for
+#'   1 - S(t), `"RMST"` for restricted mean survival time, `"incidence"` for the
+#'   person-time rate, or `"quantile"` for the first time S(t) falls to `q`.
+#' @param cluster Column name identifying the subject. Under cloning this is the
+#'   patient id, not the clone id.
+#' @param weights Weight column name, or `NULL` for an unweighted analysis.
+#' @param outcome Column name for the 0/1 outcome indicator, 1 being the event.
+#' @param time_start Column name for interval start time.
+#' @param time_stop Column name for interval stop time.
+#' @param arm Column name for treatment arm.
+#' @param horizon Time point for `"survival"`, `"cumhaz"` and `"risk"`.
+#' @param tau Restriction time for `"RMST"`.
+#' @param q Survival level for `"quantile"`; `0.5` gives median survival.
+#' @param per Person-time denominator for `"incidence"`.
+#' @param conf_level Coverage of the confidence interval.
+#'
+#' @returns A data frame with columns `est`, `se`, `lcl`, `ucl` and `p`, one row
+#'   per arm followed by the difference and ratio rows.
+#'
+#' @details
+#' Complements [emul_estimate()]: that function fits a model and returns the
+#' model object, this one returns the marginal effect measure itself with an
+#' interval. Both take the output of [weight_cases()] and the same
+#' column-naming arguments.
+#'
+#' The work runs in three layers, each callable on its own:
+#' [emul_effect_point()] for the estimate, [emul_effect_variance()] for the
+#' robust sandwich standard error, and [emul_confint()] for the interval and
+#' test. One variance object can therefore serve several intervals, and either
+#' layer can be replaced without touching the other.
+#'
+#' Little is validated, so three things are assumed and none of them raise
+#' an error when they are wrong, only a plausible number. `outcome` must be
+#' coded 0/1 with 1 the event, so survival's 1 = censored, 2 = event coding has
+#' to be recoded first. `cluster` must be the patient id, since a clone id makes
+#' a patient's two clones look independent and shrinks every contrast standard
+#' error. And intervals of one subject within one arm must be disjoint.
+#'
+#' The one thing that is checked is the time the measure is read at: a
+#' `horizon` or `tau` with no events before it, or past the end of follow-up,
+#' warns rather than quietly returning zero or an extrapolation.
+#'
+#' Weights are treated as known: no standard error here accounts for having
+#' estimated the weight model, which is conservative for a correctly specified
+#' IPTW. A quantile exists only where S(t) reaches `q` inside follow-up; where
+#' it does not the estimate is `NA`, and where it does but little of the risk
+#' set is left beneath it, the interval covers less than it claims.
+#'
+#' @export
+#' @examples
+#' set.seed(1)
+#' n <- 200
+#' treated <- rbinom(n, 1, 0.5)
+#' end <- pmin(rexp(n, ifelse(treated == 1, 0.10, 0.17)), runif(n, 2, 9))
+#' died <- as.integer(end < 9)
+#' rows <- pmax(ceiling(end), 1)
+#' dat <- data.frame(
+#'   id = rep(seq_len(n), rows),
+#'   Tstart = sequence(rows) - 1,
+#'   Tstop = pmin(sequence(rows), end[rep(seq_len(n), rows)]),
+#'   arms = factor(treated[rep(seq_len(n), rows)], 0:1, c("Control", "Surgery")),
+#'   weight_Cox = runif(sum(rows), 0.5, 2)
+#' )
+#' dat$outcome <- as.integer(dat$Tstop == end[dat$id] & died[dat$id] == 1)
+#'
+#' emul_effect(dat, "risk", weights = "weight_Cox", horizon = 5)
+#' emul_effect(dat, "RMST", weights = "weight_Cox", tau = 6)
+#' emul_effect(dat, "incidence", weights = "weight_Cox", per = 100)
+#'
+#' est <- emul_effect_point(dat, "risk", weights = "weight_Cox", horizon = 5)
+#' v <- emul_effect_variance(dat, "risk", weights = "weight_Cox", horizon = 5)
+#' emul_confint(est, v$se)
+emul_effect <- function(
+  data,
+  method = c("survival", "cumhaz", "risk", "RMST", "incidence", "quantile"),
+  cluster = "id",
+  weights = NULL,
+  outcome = "outcome",
+  time_start = "Tstart",
+  time_stop = "Tstop",
+  arm = "arms",
+  horizon = NULL,
+  tau = NULL,
+  q = 0.5,
+  per = 1000,
+  conf_level = 0.95
+) {
+  method <- match.arg(method)
+
+  estimate <- emul_effect_point(
+    data, method, cluster, weights, outcome, time_start, time_stop, arm,
+    horizon, tau, q, per
+  )
+  variance <- emul_effect_variance(
+    data, method, cluster, weights, outcome, time_start, time_stop, arm,
+    horizon, tau, q, per
+  )
+
+  emul_confint(estimate, variance$se, conf_level)
+}
+
+#' Point estimate of a marginal effect measure
+#'
+#' @inheritParams emul_effect
+#'
+#' @returns A named numeric vector, one entry per arm. With exactly two arms the
+#'   difference and the ratio are appended, always as arm 2 against arm 1.
+#' @export
+emul_effect_point <- function(
+  data,
+  method = c("survival", "cumhaz", "risk", "RMST", "incidence", "quantile"),
+  cluster = "id",
+  weights = NULL,
+  outcome = "outcome",
+  time_start = "Tstart",
+  time_stop = "Tstop",
+  arm = "arms",
+  horizon = NULL,
+  tau = NULL,
+  q = 0.5,
+  per = 1000
+) {
+  method <- match.arg(method)
+  fitted <- arm_curves(
+    data, cluster, weights, outcome, time_start, time_stop, arm
+  )
+  curves <- fitted$curves
+  check_time(fitted$data, if (identical(method, "RMST")) tau else horizon)
+
+  value <- switch(
+    method,
+    survival = vapply(curves, curve_value, numeric(1), "surv", horizon),
+    cumhaz = vapply(curves, curve_value, numeric(1), "cumhaz", horizon),
+    risk = 1 - vapply(curves, curve_value, numeric(1), "surv", horizon),
+    RMST = vapply(curves, restricted_mean, numeric(1), tau),
+    incidence = vapply(fitted$parts, function(part) {
+      per * sum(part$weight * part$event) /
+        sum(part$weight * (part$tstop - part$tstart))
+    }, numeric(1)),
+    quantile = vapply(curves, curve_quantile, numeric(1), q)
+  )
+
+  label_estimates(value, method, names(curves), q)
+}
+
+#' Robust sandwich variance of a marginal effect measure
+#'
+#' Standard errors in closed form from subject-level influence functions,
+#' `Var(theta) = sum_i psi_i^2`. Contrasts are exact rather than assumed
+#' independent, because the influence contributions of both arms live on the
+#' same subject.
+#'
+#' @inheritParams emul_effect
+#'
+#' @returns A list with `se`, whose names match [emul_effect_point()], and
+#'   `influence`, the subjects by estimands matrix of influence contributions.
+#'   A custom contrast is column arithmetic on `influence`.
+#' @export
+emul_effect_variance <- function(
+  data,
+  method = c("survival", "cumhaz", "risk", "RMST", "incidence", "quantile"),
+  cluster = "id",
+  weights = NULL,
+  outcome = "outcome",
+  time_start = "Tstart",
+  time_stop = "Tstop",
+  arm = "arms",
+  horizon = NULL,
+  tau = NULL,
+  q = 0.5,
+  per = 1000
+) {
+  method <- match.arg(method)
+  fitted <- arm_curves(
+    data, cluster, weights, outcome, time_start, time_stop, arm
+  )
+  parts <- fitted$parts
+  curves <- fitted$curves
+  ids <- fitted$ids
+
+  hazard <- function(at) {
+    vapply(parts, function(part) {
+      influence_contributions(part, ids, part$event == 1, at)
+    }, numeric(length(ids)))
+  }
+
+  block <- switch(
+    method,
+    survival = {
+      surv <- vapply(curves, curve_value, numeric(1), "surv", horizon)
+      list(psi = sweep(hazard(horizon), 2L, -surv, "*"), value = surv)
+    },
+    cumhaz = list(
+      psi = hazard(horizon),
+      value = vapply(curves, curve_value, numeric(1), "cumhaz", horizon)
+    ),
+    risk = {
+      surv <- vapply(curves, curve_value, numeric(1), "surv", horizon)
+      list(psi = sweep(hazard(horizon), 2L, surv, "*"), value = 1 - surv)
+    },
+    RMST = {
+      rmst <- vapply(curves, restricted_mean, numeric(1), tau)
+      list(
+        psi = vapply(seq_along(parts), function(i) {
+          -influence_contributions(
+            parts[[i]], ids, parts[[i]]$event == 1, tau,
+            function(times) rmst[i] - restricted_mean(curves[[i]], times)
+          )
+        }, numeric(length(ids))),
+        value = rmst
+      )
+    },
+    incidence = list(
+      psi = vapply(parts, function(part) {
+        events <- rowsum(part$weight * part$event, part$id, reorder = FALSE)
+        exposure <- rowsum(
+          part$weight * (part$tstop - part$tstart), part$id, reorder = FALSE
+        )
+        total <- sum(exposure)
+        contributions <- stats::setNames(numeric(length(ids)), ids)
+        contributions[rownames(events)] <- per *
+          (events[, 1L] - sum(events) / total * exposure[, 1L]) / total
+        contributions
+      }, numeric(length(ids))),
+      value = vapply(parts, function(part) {
+        per * sum(part$weight * part$event) /
+          sum(part$weight * (part$tstop - part$tstart))
+      }, numeric(1))
+    ),
+    quantile = {
+      theta <- vapply(curves, curve_quantile, numeric(1), q)
+      list(
+        psi = vapply(seq_along(parts), function(i) {
+          part <- parts[[i]]
+          density <- curve_density(
+            curves[[i]], theta[i], part$tstop[part$event == 1]
+          )
+          -curve_value(curves[[i]], "surv", theta[i]) *
+            influence_contributions(part, ids, part$event == 1, theta[i]) /
+            density
+        }, numeric(length(ids))),
+        value = theta
+      )
+    }
+  )
+
+  psi <- contrast_influence(block$psi, block$value, method)
+  colnames(psi) <- names(
+    label_estimates(block$value, method, names(curves), q)
+  )
+
+  list(se = sqrt(colSums(psi^2)), influence = psi)
+}
+
+#' Confidence interval and Wald test for a marginal effect measure
+#'
+#' Ratios are put on the log scale and tested against 1, differences are tested
+#' against 0. A per-arm absolute quantity has no meaningful null value, so its
+#' p is `NA` rather than a test against zero.
+#'
+#' @param estimate Named numeric vector from [emul_effect_point()].
+#' @param se Named numeric vector of standard errors, the `se` element of
+#'   [emul_effect_variance()].
+#' @param conf_level Coverage of the confidence interval.
+#'
+#' @returns A data frame with columns `est`, `se`, `lcl`, `ucl` and `p`.
+#' @export
+emul_confint <- function(estimate, se, conf_level = 0.95) {
+  terms <- names(estimate)
+  se <- se[terms]
+  contrast <- contrast_names()
+  ratio <- terms %in% contrast["ratio", ]
+  difference <- terms %in% contrast["difference", ]
+  z <- stats::qnorm(1 - (1 - conf_level) / 2)
+
+  lcl <- estimate - z * se
+  ucl <- estimate + z * se
+  lcl[ratio] <- exp(log(estimate[ratio]) - z * se[ratio] / estimate[ratio])
+  ucl[ratio] <- exp(log(estimate[ratio]) + z * se[ratio] / estimate[ratio])
+
+  null <- rep(NA_real_, length(terms))
+  null[difference] <- 0
+  null[ratio] <- 1
+  statistic <- (estimate - null) / se
+  statistic[ratio] <- log(estimate[ratio]) / (se[ratio] / estimate[ratio])
+
+  data.frame(
+    est = unname(estimate),
+    se = unname(se),
+    lcl = unname(lcl),
+    ucl = unname(ucl),
+    p = unname(2 * stats::pnorm(-abs(statistic))),
+    row.names = terms
+  )
+}