A Case Study Using the Beta-Danish Distribution
Source:vignettes/betadanish-case-study.Rmd
betadanish-case-study.RmdIntroduction
This vignette demonstrates a typical survival analysis workflow using the BetaDanish package.
library(BetaDanish)
#> BetaDanish 0.3.0: see ?BetaDanish for help.
library(survival)
#>
#> 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 1Fitting 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.990, 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.897875 -1.073246 2.446423 0.7647 0.444461
#> b 4.078124 1.491051 1.155665 7.000584 2.7351 0.006237 **
#> c 2.196484 2.901329 -3.490120 7.883088 0.7571 0.449013
#> k 0.082975 0.084062 -0.081788 0.247737 0.9871 0.323613
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> ---
#> Log-Likelihood: -409.9137
#> AIC: 827.8274 | BIC: 839.2356Three-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
