---
title: "Using your own data generator"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Using your own data generator}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>"
)
```

## Overview

`simdata_fast()` generates piecewise-exponential survival and dropout times,
subgroups, multi-arm trials, and correlated endpoints from an illness-death
model. Other data-generating models, such as Weibull or cure-model survival,
copula-based dependence, or resampling from historical data, are easy to write
in a few lines of R. The analysis side of FastSurvival does not depend on how
the data were generated: `cutoff_fast()`, `analysis_fast()`, `switch_fast()`,
and `simsummary_fast()` read a data frame with a fixed set of columns. This
vignette generates trials from a Weibull model with a cure fraction in the
treatment group, a model that `simdata_fast()` does not provide, and runs the
same workflow as for data from `simdata_fast()`.

The columns that the functions read are

| Function | Required columns |
|----------|------------------|
| `cutoff_fast()` | `sim`, `accrual_time`, and the time and event columns named by `tte.col` and `event.col` (`tte` and `event` by default) |
| `analysis_fast()`, `pairwise_fast()` | `sim`, `group`, `accrual_time`, `tte`, `event` |
| `switch_fast()` | `sim`, `group`, `accrual_time`, and the latent times `surv_time` and `dropout_time` |

where `tte` is the observed time from accrual, `event` is 1 for an event and 0
for censoring, and `accrual_time` is the calendar time of enrollment. The
latent columns are needed only by `switch_fast()`, which changes the survival
time after the switch and recomputes the observed columns.

```{r load}
library(FastSurvival)
```

## A Weibull model with a cure fraction

The control group has Weibull survival with shape 1.3 and scale 16 months. In
the treatment group 15% of the subjects are cured and the others have Weibull
survival with the same shape and scale 18 months, so the survival curves
separate slowly and the treatment curve levels off at 0.15. Subjects are
enrolled uniformly over 12 months and drop out at a constant rate of 0.01 per
month. The generator below writes the columns listed above for all simulated
trials at once.

```{r generator}
gen_cure_weibull <- function(nsim, n_per, accrual, shape, scale, cure,
                             drop_rate, seed) {
  set.seed(seed)
  n_tot <- 2L * n_per
  N     <- nsim * n_tot
  sim   <- rep(seq_len(nsim), each = n_tot)
  group <- rep(rep(1:2, each = n_per), times = nsim)
  accrual_time <- stats::runif(N, 0, accrual)
  surv_time    <- stats::rweibull(N, shape = shape, scale = scale[group])
  surv_time[stats::runif(N) < cure[group]] <- Inf   # cured subjects
  dropout_time <- stats::rexp(N, rate = drop_rate)
  tte <- pmin(surv_time, dropout_time)
  data.frame(sim, group, accrual_time, surv_time, dropout_time, tte,
             event = as.integer(surv_time <= dropout_time),
             calendar_time = accrual_time + tte)
}

shape <- 1.3
scale <- c(16, 18)
cure  <- c(0, 0.15)

dat <- gen_cure_weibull(nsim = 2000, n_per = 200, accrual = 12,
                        shape = shape, scale = scale, cure = cure,
                        drop_rate = 0.01, seed = 2026)
head(dat)
```

A quick check of the generator compares the proportion of latent survival
times beyond a few time points with the survival function of the model,
`cure + (1 - cure) * exp(-(t / scale)^shape)`.

```{r check-generator}
surv_model <- function(t, g) cure[g] + (1 - cure[g]) * exp(-(t / scale[g])^shape)
t_chk <- c(6, 12, 24, 36)
data.frame(
  t               = t_chk,
  control_sim     = sapply(t_chk, function(t) mean(dat$surv_time[dat$group == 1] > t)),
  control_model   = surv_model(t_chk, 1),
  treatment_sim   = sapply(t_chk, function(t) mean(dat$surv_time[dat$group == 2] > t)),
  treatment_model = surv_model(t_chk, 2)
)
```

## Analysis times

The interim analysis takes place at 150 events or at month 24, whichever comes
first, and the final analysis at 280 events, but not earlier than 9 months
after the interim and not later than month 48. `cutoff_fast()` returns the
calendar time of both looks for every simulated trial.

```{r cutoffs}
cut <- cutoff_fast(dat, event.looks = c(150, 280), max.time = c(24, 48),
                   min.gap = c(NA, 9))
head(cut)
colMeans(cut)
```

## Log-rank and max-combo tests

Because the treatment effect builds up late and ends in a plateau, the
max-combo test of Fleming-Harrington weights is compared with the log-rank
test. The nominal one-sided levels 0.002 and 0.024 are close to those of an
O'Brien-Fleming-type spending function at these looks and are used for
illustration. With `mc.alpha` set to these levels, the max-combo p-value is
integrated only when the Bonferroni bounds do not decide the comparison with
the level, which leaves the decisions unchanged.

```{r analysis}
alpha_look <- c(0.002, 0.024)
set.seed(1)
res <- analysis_fast(dat, control = 1, cutoff.looks = cut,
                     stat = c("logrank", "maxcombo"), side = 1,
                     mc.alpha = alpha_look)

oc_lr <- simsummary_fast(res, p.col = "logrank.p",  alpha = alpha_look)
oc_mc <- simsummary_fast(res, p.col = "maxcombo.p", alpha = alpha_look)
data.frame(
  test             = c("Log-rank", "Max-combo"),
  reject_interim   = c(oc_lr[oc_lr$look == "1", "prob.stop.efficacy"],
                       oc_mc[oc_mc$look == "1", "prob.stop.efficacy"]),
  power            = c(oc_lr[oc_lr$look == "overall", "cum.reject"],
                       oc_mc[oc_mc$look == "overall", "cum.reject"])
)
mean(res$maxcombo.p.exact, na.rm = TRUE)
```

The last value is the share of the max-combo p-values that required the
multivariate normal integral.

## Crossover after the interim analysis

`switch_fast()` works on the same data. Here every control subject who is still
on study at the interim analysis switches to the experimental treatment, which
multiplies the remaining survival time by 1.5. The final analysis is then
triggered by the events of the modified data. The outcomes observed by the
interim analysis are unchanged, so the interim cutoffs and statistics are the
same as before.

```{r switching}
sw <- switch_fast(dat, group = 1, when = "cutoff",
                  cutoff = cut[, 1, drop = FALSE], aft.factor = 1.5)
cut_sw <- cutoff_fast(sw, event.looks = c(150, 280), max.time = c(24, 48),
                      min.gap = c(NA, 9))
all.equal(cut_sw[, 1], cut[, 1])

res_sw <- analysis_fast(sw, control = 1, cutoff.looks = cut_sw,
                        stat = "logrank", side = 1)
all.equal(res_sw$logrank.z[res_sw$look == 1], res$logrank.z[res$look == 1])

oc_sw <- simsummary_fast(res_sw, p.col = "logrank.p", alpha = alpha_look)
c(log_rank_power_without_switching = oc_lr[oc_lr$look == "overall", "cum.reject"],
  log_rank_power_with_crossover    = oc_sw[oc_sw$look == "overall", "cum.reject"])
```

## Remarks

Any generator can be used in the same way, provided that it writes one row per
subject with the columns above and that the simulation identifiers group the
rows of each trial. The data do not need to be sorted, but generating them
grouped by `sim`, as `simdata_fast()` does, avoids a sorting step in the
analysis functions. For an illness-death model, `switch_fast()` reads the
columns `e1_surv_time`, `e2_surv_time`, `dropout_time`, and `intermediate`
instead, as described in its help page.