A Case Study Using the Beta-Danish Distribution

Introduction

This vignette demonstrates a typical survival analysis workflow using the BetaDanish package.

library(BetaDanish)
library(survival)
#> Warning: package 'survival' was built under R version 4.5.3
#> 
#> Attaching package: 'survival'
#> The following objects are masked from 'package:BetaDanish':
#> 
#>     leukemia, transplant
data('remission', package = 'BetaDanish')
head(remission)
#>   time status
#> 1 0.08      1
#> 2 2.09      1
#> 3 3.48      1
#> 4 4.87      1
#> 5 6.94      1
#> 6 8.66      1

Fitting the Beta-Danish model

fit <- fit_betadanish(Surv(time, status) ~ 1, data = remission, n_starts = 1)
#> Warning: The times look recorded on a grid of 0.01. The point-density
#> likelihood treats them as exact, which understates the standard errors.
#> Consider grouped = TRUE.
#> Warning: The fitted correlation between a-hat and c-hat is -0.992, so
#> effectively only the product c*a is identified. Individual estimates of a and c
#> should not be interpreted.
#> Warning: 1 starting point(s) reached a degenerate ridge and were discarded. The
#> reported fit is the best admissible optimum. If this is most of the grid, the
#> four-parameter model is a poor choice for these data.
summary(fit)
#> 
#> Call:
#> fit_betadanish(formula = Surv(time, status) ~ 1, data = remission, 
#>     n_starts = 1)
#> 
#> Beta-Danish Distribution Fit
#> Model: Full 4-Parameter Model 
#> 
#>    Estimate Std. Error Lower 95% Upper 95% z value Pr(>|z|)   
#> a  0.686589   0.980513 -1.235217  2.608395  0.7002 0.483781   
#> b  4.078124   1.529590  1.080128  7.076120  2.6662 0.007672 **
#> c  2.196484   3.184261 -4.044668  8.437637  0.6898 0.490324   
#> k  0.082975   0.092057 -0.097456  0.263405  0.9013 0.367405   
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> ---
#> Log-Likelihood: -409.9137 
#> AIC: 827.8274  | BIC: 839.2356

Three-parameter submodel

fit_sub <- fit_betadanish(Surv(time, status) ~ 1, data = remission, submodel = TRUE, n_starts = 1)
#> Warning: The times look recorded on a grid of 0.01. The point-density
#> likelihood treats them as exact, which understates the standard errors.
#> Consider grouped = TRUE.
#> Warning: b-hat is only 1.90 standard errors from 1, close to the b = 1
#> non-identifiability ridge. Consider the ED submodel (submodel = TRUE). See the
#> Identifiability section of ?fit_betadanish.
compare_models(fit, fit_sub)
#> Likelihood Ratio Test (a = 1 vs a != 1)
#> 
#>                  Model    LogLik      Chisq Df Pr(>Chisq)
#> 1   Submodel (3-param) -409.9541         NA NA         NA
#> 2 Full Model (4-param) -409.9137 0.08081129  1   0.776201

Diagnostic plots

plot(fit, type = 'survival')

plot(fit, type = 'hazard')

Interpretation

The fitted model can be used to estimate survival probabilities, hazard behavior, and overall model fit. Users should compare the Beta-Danish model with alternative lifetime distributions and inspect diagnostic plots before drawing final conclusions.