Competing Risks with the Beta-Danish Distribution

Simulated two-cause data

library(BetaDanish)
set.seed(2026)
n  <- 400
T1 <- rbetadanish(n, a = 1.2, b = 1.5, c = 1.0, k = 0.4)
T2 <- rbetadanish(n, a = 1.0, b = 2.0, c = 1.0, k = 0.2)
C  <- stats::rexp(n, 0.05)
time  <- pmin(T1, T2, C)
cause <- ifelse(time == C, 0L, ifelse(T1 <= T2, 1L, 2L))
table(cause)
#> cause
#>   0   1   2 
#>  23 198 179

Fitting the model

fit <- fit_bd_competing(time = time, cause = cause)
print(fit)
#> 
#> Call:
#> fit_bd_competing(time = time, cause = cause)
#> 
#> Beta-Danish Competing Risks Model
#> Log-Likelihood: -700.0237 
#> 
#> Cause-specific estimates:
#>              a        b        c      k
#> Cause_1 0.0054 132.5037 271.2345 1.7050
#> Cause_2 0.0630   2.4563  17.2279 0.4662

CIF comparison: Beta-Danish vs Aalen-Johansen

res <- cif_compare(fit, plot = TRUE)

Gray’s test

if (!is.null(res$gray_test)) {
  print(res$gray_test)
} else {
  cat("Gray's test was not produced.\n")
}
#> Gray's test was not produced.

Covariates (new in 0.3.0)

Each cause can carry its own coefficient vector, acting on the time scale:

\[T_j = T_{0j}\exp(x'\gamma_j)\]

so a positive coefficient lengthens time to failure from that cause. The same covariate may therefore accelerate one cause and retard another, which is the point of fitting them separately.

# The covariate column only exists when 'gammas' is supplied
d <- simulate_bd_competing_data(400, gammas = c(0.8, -0.8), seed = 1)

fit <- fit_bd_competing(d$time, d$cause, covariates = ~ x, data = d,
                        submodel = TRUE)
fit$coefficients

The intercept is dropped from the design matrix, because each cause already has its own scale parameter k that an intercept would be confounded with.

cif_betadanish() takes the covariate values at which to evaluate the cumulative incidence, defaulting to the reference subject:

cif_betadanish(fit, tvec = c(5, 20, 50), cause_idx = 1, x = 0)
cif_betadanish(fit, tvec = c(5, 20, 50), cause_idx = 1, x = 1)

A covariate effect folds exactly into the scale: scaling time by \(\lambda\) maps \(k\) to \(k/\lambda\), so no re-integration is needed.

A caution worth repeating

Independence of the latent failure times is an identifying assumption, not a testable one. Under positive dependence the working independence model biases the cause-specific cumulative incidence functions downward. Overall survival is more robust than the individual marginals; where conclusions rest on absolute CIFs, add a copula-based sensitivity analysis.