This reference covers two time-to-event / latent-process simulators. simulate_survival() draws right-censored survival data from a Cox proportional-hazards generative model with known log hazard ratios, and simulate_hmm() produces observed symbol sequences from a discrete hidden Markov model with known transition and emission matrices. Both return a saqr_sim object whose $params hold the ground truth.

simulate_survival()

Generate right-censored survival data from a Cox proportional-hazards model with a parametric baseline hazard (Weibull, exponential, or Gompertz). The data are built so that coxph(Surv(time, status) ~ ., data) recovers the true betas at large n.

Signature

args(simulate_survival)
#> function (n = 300, n_covariates = 2, betas = NULL, baseline = c("weibull", 
#>     "exponential", "gompertz"), lambda = 0.1, shape = 1, censoring_rate = 0.3, 
#>     covariate_type = c("normal", "binary"), seed = NULL) 
#> NULL

Example

The data frame has time, a status indicator (1 = event, 0 = censored), and the covariates X1 ... Xp.

r <- simulate_survival(n = 400, n_covariates = 2, seed = 1)
r
#> saqr_sim [survival]  400 x 4  (seed=1)
#>   params: betas, baseline, lambda, shape, covariate_type, realized_censoring_rate 
#>   cols:   time, status, X1, X2
head(as.data.frame(r))
r$params$betas
#>   X1   X2 
#>  0.5 -0.3
r$params$realized_censoring_rate
#> [1] 0.2975

A second configuration with binary covariates and an exponential baseline:

r2 <- simulate_survival(n = 500, betas = c(0.8, -0.5),
                        baseline = "exponential",
                        covariate_type = "binary", seed = 42)
summary(r2)
#> Simulation type: survival
#> Seed: 42
#> 
#> --- Data ---
#>       time               status            X1             X2       
#>  Min.   : 0.002791   Min.   :0.000   Min.   :0.00   Min.   :0.000  
#>  1st Qu.: 1.863159   1st Qu.:0.000   1st Qu.:0.00   1st Qu.:0.000  
#>  Median : 4.211101   Median :1.000   Median :0.00   Median :0.000  
#>  Mean   : 6.916516   Mean   :0.674   Mean   :0.47   Mean   :0.476  
#>  3rd Qu.: 9.141489   3rd Qu.:1.000   3rd Qu.:1.00   3rd Qu.:1.000  
#>  Max.   :51.034301   Max.   :1.000   Max.   :1.00   Max.   :1.000  
#> 
#> --- Parameters ---
#> List of 6
#>  $ betas                  : Named num [1:2] 0.8 -0.5
#>   ..- attr(*, "names")= chr [1:2] "X1" "X2"
#>  $ baseline               : chr "exponential"
#>  $ lambda                 : num 0.1
#>  $ shape                  : num 1
#>  $ covariate_type         : chr "binary"
#>  $ realized_censoring_rate: num 0.326

Visualising event vs. censoring times

Overlaid histograms of observed times, split by whether the observation was an event or was censored, show the censoring structure.

surv_df <- as.data.frame(r)
ev   <- surv_df$time[surv_df$status == 1]
cens <- surv_df$time[surv_df$status == 0]
brk  <- seq(0, max(surv_df$time) + 1, length.out = 25)

hist(ev, breaks = brk, col = "#2c7fb880", border = "white",
     xlab = "Observed time", main = "Event vs. censoring times",
     ylim = c(0, max(hist(ev, breaks = brk, plot = FALSE)$counts,
                     hist(cens, breaks = brk, plot = FALSE)$counts)))
hist(cens, breaks = brk, col = "#d95f0e80", border = "white", add = TRUE)
legend("topright", fill = c("#2c7fb880", "#d95f0e80"),
       legend = c("event (status = 1)", "censored (status = 0)"), bty = "n")

Optional: Kaplan–Meier and recovery with survival

When the survival package is installed, a Kaplan–Meier curve and a Cox fit confirm the generating effects. Skipped if survival is unavailable.

km <- survival::survfit(survival::Surv(time, status) ~ 1, data = surv_df)
plot(km, col = "#2c7fb8", lwd = 2, xlab = "Time", ylab = "Survival probability",
     main = "Kaplan–Meier estimate")


cox <- survival::coxph(survival::Surv(time, status) ~ X1 + X2, data = surv_df)
data.frame(true = r$params$betas, estimated = unname(coef(cox)))

simulate_hmm()

Generate observed symbol sequences from a discrete hidden Markov model. For each sequence a latent state path is drawn from the Markov chain (init, trans), and each hidden state emits an observed symbol via the emission matrix. The true latent paths are returned in $params$hidden_paths.

Signature

args(simulate_hmm)
#> function (n_sequences = 50, seq_length = 30, n_states = 2, n_symbols = 3, 
#>     trans = NULL, emission = NULL, init = NULL, seed = NULL) 
#> NULL

Example

The data are long-format: one row per (sequence_id, time) with the observed symbol. The true transition and emission matrices live in $params.

h <- simulate_hmm(n_sequences = 100, seq_length = 40, seed = 1)
h
#> saqr_sim [hmm]  4000 x 3  (seed=1)
#>   params: trans, emission, init, n_states, n_symbols, hidden_paths 
#>   cols:   sequence_id, time, symbol
head(as.data.frame(h))
h$params$trans
#>      [,1] [,2]
#> [1,]  0.8  0.2
#> [2,]  0.2  0.8
h$params$emission
#>      [,1] [,2] [,3]
#> [1,]  0.8  0.1  0.1
#> [2,]  0.1  0.8  0.1

A second, explicit three-state model:

tr <- matrix(c(0.8, 0.1, 0.1,
               0.1, 0.8, 0.1,
               0.1, 0.1, 0.8), nrow = 3, byrow = TRUE)
em <- matrix(c(0.7, 0.2, 0.1,
               0.1, 0.7, 0.2,
               0.2, 0.1, 0.7), nrow = 3, byrow = TRUE)
h2 <- simulate_hmm(n_sequences = 80, seq_length = 50, n_states = 3,
                   n_symbols = 3, trans = tr, emission = em, seed = 42)
summary(h2)
#> Simulation type: hmm
#> Seed: 42
#> 
#> --- Data ---
#>   sequence_id         time          symbol     
#>  Min.   : 1.00   Min.   : 1.0   Min.   :1.000  
#>  1st Qu.:20.75   1st Qu.:13.0   1st Qu.:1.000  
#>  Median :40.50   Median :25.5   Median :2.000  
#>  Mean   :40.50   Mean   :25.5   Mean   :2.002  
#>  3rd Qu.:60.25   3rd Qu.:38.0   3rd Qu.:3.000  
#>  Max.   :80.00   Max.   :50.0   Max.   :3.000  
#> 
#> --- Parameters ---
#> List of 6
#>  $ trans       : num [1:3, 1:3] 0.8 0.1 0.1 0.1 0.8 0.1 0.1 0.1 0.8
#>  $ emission    : num [1:3, 1:3] 0.7 0.1 0.2 0.2 0.7 0.1 0.1 0.2 0.7
#>  $ init        : num [1:3] 0.333 0.333 0.333
#>  $ n_states    : int 3
#>  $ n_symbols   : int 3
#>  $ hidden_paths: int [1:80, 1:50] 1 3 3 1 1 3 3 3 2 3 ...
dim(h2$params$hidden_paths)
#> [1] 80 50
h2$params$hidden_paths[1:3, 1:10]
#>      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,]    1    3    3    3    3    3    3    2    2     2
#> [2,]    3    3    3    3    1    1    1    1    1     1
#> [3,]    3    3    3    3    3    3    3    3    3     1