Introduction to the BetaDanish Package

1. Introduction

The BetaDanish package provides a comprehensive suite of tools for survival and reliability analysis using the Beta-Danish distribution.

Classical lifetime models, such as the Weibull, Gamma, or Log-Normal distributions, are often limited in their hazard shape flexibility. They typically force the hazard rate to be strictly monotonic or, at best, unimodal. The Beta-Danish distribution overcomes these limitations by accommodating decreasing, increasing, unimodal (upside-down bathtub), and bathtub-shaped hazard rates within a single, unified four-parameter family.

This package implements the core mathematical functions, robust Maximum Likelihood Estimation (MLE) for complete and right-censored data, and advanced modules for Accelerated Failure Time (AFT) regression, cure-rate models, and competing risks.

2. Theoretical Background

The Beta-Danish distribution is constructed by applying the beta-generated transformation (Eugene et al., 2002) to a baseline distribution. The baseline used here is the two-parameter exponentiated log-logistic distribution.

Let \(a, b, c, k > 0\) be the parameters of the distribution. * \(k\) is the scale parameter. * \(c\) is the baseline shape parameter. * \(a\) and \(b\) are the beta-generator shape parameters controlling skewness and tail weight.

The Cumulative Distribution Function (CDF) is given by the regularized incomplete beta function: \[ F(t) = I_{G(t)}(a, b) \] where \(G(t) = \left( \frac{k t}{1 + k t} \right)^c\).

The Probability Density Function (PDF) is: \[ f(t) = \frac{c k}{B(a, b)} \frac{(k t)^{ca - 1}}{(1 + k t)^{ca + 1}} \left[ 1 - \left( \frac{k t}{1 + k t} \right)^c \right]^{b - 1} \]

The 3-Parameter Submodel

By fixing \(a = 1\), the distribution reduces to a highly tractable 3-parameter submodel (the Complementary Exponentiated Danish distribution), which is particularly useful for regression and cure-rate modeling to ensure parameter identifiability.

3. Core Distribution Functions

The package provides the standard d/p/q/r/h functions. All internal calculations are performed in log-space to ensure numerical stability, especially in the tails.

library(BetaDanish)

# Calculate the hazard rate at time t = 2
hbetadanish(x = 2, a = 1.5, b = 2.0, c = 3.0, k = 0.5)

# Generate 100 random survival times
set.seed(2026)
sim_data <- rbetadanish(n = 100, a = 1.5, b = 2.0, c = 3.0, k = 0.5)

4. Fitting Models to Data

The fit_betadanish() function is the core MLE engine. It uses a multi-start optimization strategy to ensure global convergence and strictly enforces positivity constraints via log-scale reparameterization.

Example 1: Uncensored Data (Bladder Cancer Remission)

The remission dataset contains the remission times of 128 bladder cancer patients. This data exhibits a unimodal/decreasing hazard rate.

# Load the built-in dataset
data("remission", package = "BetaDanish")

# Fit the full 4-parameter model
fit_full <- fit_betadanish(survival::Surv(time, status) ~ 1, data = remission)
summary(fit_full)

# Fit the 3-parameter submodel (a = 1)
fit_sub <- fit_betadanish(survival::Surv(time, status) ~ 1, data = remission, submodel = TRUE)

# Compare the nested models using a Likelihood Ratio Test (LRT)
compare_models(fit_full, fit_sub)

You can easily generate publication-quality diagnostic plots:

# Generates Survival, Hazard, Density, P-P, and Q-Q plots
plot(fit_full, type = "all")

Example 2: Right-Censored Data (Bone Marrow Transplant)

The transplant dataset contains survival times for 91 patients, including right-censored observations. The fit_betadanish() function natively handles survival::Surv objects.

data("transplant", package = "BetaDanish")

# Fit the model to censored data
fit_censored <- fit_betadanish(survival::Surv(time, status) ~ 1, data = transplant)
summary(fit_censored)

# Plot the Kaplan-Meier curve overlaid with the Beta-Danish fit
plot(fit_censored, type = "survival")

5. Advanced Modeling: Cure Rates

In many clinical datasets (like the transplant data), a proportion of patients may be “cured” and will never experience the event. Standard survival models force the survival curve to zero, which misrepresents the data.

The BetaDanish package provides the fit_bd_cure() function to fit both Mixture and Promotion-Time (Non-Mixture) cure models.

# Fit a mixture cure model
# Latency (time to event for susceptible) has no covariates (~ 1)
# Incidence (probability of being susceptible) depends on treatment group (~ group)
cure_fit <- fit_bd_cure(
  formula_aft = survival::Surv(time, status) ~ 1,
  formula_cure = ~ group,
  data = transplant,
  type = "mixture"
)

summary(cure_fit)

6. Automated Benchmarking

To justify the use of the Beta-Danish distribution, you can benchmark it against classical distributions (Weibull, Gamma, Log-Normal, etc.) using the compare_distributions() function. (Note: This requires the flexsurv package).

# Returns a ranked table of AIC, BIC, and Log-Likelihoods
compare_distributions(fit_full)

7. Conclusion

The BetaDanish package offers a robust, flexible, and user-friendly environment for advanced survival analysis. By unifying standard MLE fitting, comprehensive diagnostics, and advanced modules (AFT, Cure, Competing Risks) under a single framework, it serves as a powerful tool for statisticians and applied researchers alike.

New in 0.3.0

Three groups of functions were added after the 0.2.0 release.

Structural properties

p <- list(a = 1.5, b = 5, c = 2, k = 1)

bd_moment_summary(p$a, p$b, p$c, p$k)
#>       mean   variance         sd   skewness   kurtosis 
#>  1.0517578  0.7961493  0.8922720  3.6708237 49.3171010 
#> attr(,"condition")
#> [1] "E(Z^r) is finite if and only if b > r"
bd_entropy_shannon(p$a, p$b, p$c, p$k)
#> [1] 0.9104078
bd_stress_strength(strength = p, stress = p)   # identical laws: one half
#> [1] 0.5

Moments exist only when b > r, and the package enforces it rather than returning a plausible-looking number from a truncated sum:

bd_moments(1:4, a = 1.5, b = 3, c = 2, k = 1)   # the fourth is Inf
#> [1] 1.84375 8.87500     Inf     Inf
bd_tail_index(a = 1.5, b = 3, c = 2, k = 1)$moment_condition
#> [1] "E(Z^r) finite if and only if b > r"

Shape, before and after fitting

bd_ttt_plot() reads the hazard shape from the data alone, with no model assumed. bd_hazard_shape() does the same from a fitted parameter set, so the two can be compared.

data(guinea_pig)
ttt <- bd_ttt_plot(guinea_pig$time)

attr(ttt, "shape")
#> [1] "unimodal (upside-down bathtub)"

Interval estimation

All four parameters are positive, so a symmetric Wald interval on the natural scale can cross zero. bd_wald_ci() builds on the log scale instead, and bd_profile_ci() inverts the likelihood ratio test. When the profile never returns below the critical threshold there is no finite upper bound, and the honest report is a lower bound rather than a point estimate with an interval around it.

fit <- fit_betadanish(survival::Surv(time, status) ~ 1, data = remission,
                      submodel = TRUE)
bd_wald_ci(fit)
bd_profile_plot(bd_profile_ci(fit, "b"))

# For the four-parameter model, report the identified composite
bd_identified_coef(fit_betadanish(survival::Surv(time, status) ~ 1,
                                  data = remission))

Simulation studies

bd_simulation_study(), bd_simulation_cure() and bd_simulation_competing() report bias, RMSE, the ratio of mean estimated standard error to empirical standard deviation, and Wald coverage. Read the last two together: a ratio well below one means the asymptotic standard errors are optimistic at that sample size, and coverage falls short as a consequence.

bd_simulation_study(n = c(50, 100, 200), n_sim = 500,
                    truth = c(b = 3, c = 2, k = 0.5), submodel = TRUE)