Fits the Beta-Danish distribution by maximum likelihood. Complete and right-censored samples are both supported, via a `survival::Surv` response.
Usage
fit_betadanish(
formula,
data,
submodel = FALSE,
n_starts = 10,
method = "BFGS",
check_identifiability = TRUE,
penalty = 0,
penalty_center = NULL,
grouped = FALSE,
delta = NULL
)Arguments
- formula
A formula whose left-hand side is a `Surv` object. Use `~ 1` for a model without covariates.
- data
A data frame containing the variables in `formula`.
- submodel
Logical; if `TRUE`, fits the three-parameter Exponentiated Danish (ED) submodel by fixing `a = 1`.
- n_starts
Integer; number of random starting points for the multi-start optimisation. Default 10.
- method
Character; optimisation method passed to `maxLik::maxLik`.
- check_identifiability
Logical; if `TRUE` (default), issue warnings when the fit lands in a region where the parameters are weakly identified.
- penalty
Non-negative ridge penalty applied on the log-parameter scale. `0` (default) is ordinary maximum likelihood. A small positive value stabilises the fit when the likelihood is nearly flat, at the cost of bias toward `penalty_center`.
- penalty_center
Numeric vector on the log-parameter scale toward which the penalty shrinks, or `NULL` (default) to shrink toward the best unpenalised start.
- grouped
Logical; if `TRUE`, use the grouped (interval) likelihood appropriate to times recorded on a coarse grid. Default `FALSE`.
- delta
Recording increment for `grouped = TRUE`. `NULL` (default) infers it from the spacing of the observed times.
Value
An object of S3 class `"betadanish"` with components including `coefficients`, `logLik`, `vcov`, `npar`, `nobs`, `convergence` and `diagnostics`.
Details
Optimisation is carried out on log-transformed parameters so that positivity is enforced without constraints; estimates and the variance-covariance matrix are returned on the natural scale, the latter via the delta method.
Grouped data
Survival times are often recorded on a grid – whole days, whole months – and the point-density likelihood is not appropriate for them. It treats a rounded value as an exact observation, which can overstate the information in the sample and make uncertainty estimates too small. With `grouped = TRUE` an event recorded at \(t\) contributes \(\log\{F(t + \delta/2) - F(t - \delta/2)\}\) instead of \(\log f(t)\); censored observations are unchanged. The cell probability falls back to \(f(t)\delta\) only where the difference of two nearly equal distribution values has cancelled to zero.
`read_survival_data()` reports an inferred `grid_step` for exactly this purpose, and this function warns when the times look grid-recorded but `grouped = FALSE`.
The grouped likelihood is a harder surface to optimise than the point-density one. On the same data and the same starting grid, the optimiser can stop where the observed information is not positive definite, and the delta-method variances then come back non-positive. Check `fit$diagnostics$vcov_singular` before trusting a standard error from a grouped fit; leaving `check_identifiability = TRUE` will warn about it automatically.
Penalised fitting
`penalty > 0` adds \(\lambda \sum (\theta - \mu)^2\) on the log-parameter scale. This is worth reaching for when the likelihood is flat along the \((a, c)\) direction and the unpenalised optimiser wanders, but it is a deliberate bias: the estimates are shrunk toward `penalty_center`.
The reported `logLik` is always the **unpenalised** log-likelihood evaluated at the penalised estimate, so that AIC, BIC and likelihood ratio tests remain comparable across fits. The objective actually maximised is stored separately as `penalised_logLik`. Treat the degrees of freedom as nominal: shrinkage reduces the effective number of parameters, so information criteria are conservative under penalisation.
Identifiability
The four-parameter model is not uniformly well identified, and a converged fit is not by itself evidence that it is. Two regions warrant care.
* **The \(b = 1\) ridge.** At \(b = 1\) the beta generator collapses and the model is non-identifiable. A fit with \(\hat b\) within about two standard errors of one lies close to that ridge; the likelihood is nearly flat along it, so the individual estimates carry little information even though the fitted survival curve may look excellent. * **Lower-tail \((a, c)\) confounding.** Near the lower tail, \(a\) and \(c\) enter almost exclusively through the product \(ca\), so the expected Fisher information is close to singular in that direction. A fitted correlation between \(\hat a\) and \(\hat c\) above about 0.95 in absolute value indicates that only the product is being estimated.
In either case the ED submodel (`submodel = TRUE`) is usually the honest report, and a likelihood ratio test via [compare_models()] will normally fail to reject it. Set `check_identifiability = FALSE` to silence the warnings once you have satisfied yourself that they are understood.
Examples
# \donttest{
set.seed(123)
sim_time <- rbetadanish(150, a = 1.5, b = 3, c = 2, k = 0.5)
sim_status <- rbinom(150, 1, 0.85)
dat <- data.frame(time = sim_time, status = sim_status)
fit <- fit_betadanish(survival::Surv(time, status) ~ 1, data = dat)
#> Warning: The fitted correlation between a-hat and c-hat is -0.989, so effectively only the product c*a is identified. Individual estimates of a and c should not be interpreted.
summary(fit)
#>
#> Call:
#> fit_betadanish(formula = survival::Surv(time, status) ~ 1, data = dat)
#>
#> Beta-Danish Distribution Fit
#> Model: Full 4-Parameter Model
#>
#> Estimate Std. Error Lower 95% Upper 95% z value Pr(>|z|)
#> a 1.27620 2.22797 -3.09063 5.64302 0.5728 0.566776
#> b 3.41176 1.08464 1.28587 5.53765 3.1455 0.001658 **
#> c 2.04265 3.54925 -4.91387 8.99917 0.5755 0.564942
#> k 0.33468 0.32404 -0.30044 0.96979 1.0328 0.301687
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> ---
#> Log-Likelihood: -302.6948
#> AIC: 613.3895 | BIC: 625.4321
fit_sub <- fit_betadanish(survival::Surv(time, status) ~ 1, data = dat,
submodel = TRUE)
compare_models(fit, fit_sub)
#> Likelihood Ratio Test (a = 1 vs a != 1)
#>
#> Model LogLik Chisq Df Pr(>Chisq)
#> 1 Submodel (3-param) -302.7004 NA NA NA
#> 2 Full Model (4-param) -302.6948 0.01119187 1 0.9157477
# }