This release brings the package into line with the final analysis code of the doctoral thesis on which it is based (Ahmad, 2026). Where a routine of the package and of the thesis estimate the same quantity, they now use the same method.
fit_betadanish() uses the two-stage estimator of the
thesis. Stage 1 fits the Exponentiated Danish (ED) submodel from four
starts with k = 1/median(t); Stage 2 fits the
four-parameter model from seven starts along the (a, c)
ridge through the Stage-1 optimum,
(f, b_ED, c_ED / f, k_ED) for f in
{1/30, 1/8, 1/3, 1, 3, 8, 30}. The start f = 1
is the ED optimum, so the four-parameter log-likelihood can never fall
below the ED one. The deterministic grid and random starts of 0.3.0 are
kept as further candidates, drawn in the same order. For complete,
ungrouped samples the analytic score (thesis Appendix B.1) is used as
the gradient. The selected optimum is then polished to a relative
tolerance of 1e-12; the default of maxLik,
about 1.5e-8, could leave the log-likelihood some
1e-5 short of the maximum. The fitted object gains
submodel_logLik.
Parameter box. The shape parameters are confined
to the thesis box [1e-4, 1e4] (new argument
box); the rate k gets the same box widened to
follow the unit of time. A fit within 2% of the edge, on the log scale,
is flagged at_boundary = TRUE: its log-likelihood is a
supremum, so AIC, BIC, the standard errors and
the goodness-of-fit statistics are NA,
summary() and report_betadanish() say so, and
a warning explains it. An interior fit with a < 0.01 or
c > 100 is flagged as weakly identified on the
(a, c) ridge but keeps its information criteria. New
components at_boundary, edge_par,
flag, box and bounds.
fit_bd_cure() gains submodel. With
submodel = FALSE it fits the four-parameter Beta-Danish
kernel (adding log_a), started from the ED cure optimum and
from the (a, c) ridge through it, and refuses a fit whose
log-likelihood is below that of its own ED submodel. Cure fits now carry
npar, nobs, AIC,
BIC, logLik_submodel and
separated, the last flagging a regression coefficient above
15 in absolute value (near-separation). summary() reports
the likelihood ratio test of a = 1. The default,
submodel = TRUE, gives the same fit as before.
fit_bd_competing() takes its starting values from
the cause-wise univariate fits only when they are admissible for the
joint model. A cause-wise fit on the edge of the parameter box, or with
a shape above 500, is replaced by that cause’s ED fit
(a = 1), so the joint search is never seeded on the
degenerate ridge.
gof_betadanish() now reports the Kolmogorov-Smirnov,
Cramer-von Mises and Anderson-Darling statistics, computed with the
Kaplan-Meier estimator on the probability-integral-transform scale.
Without censoring they are exactly the classical statistics, so complete
and censored samples are on one scale. New components
CvM_Statistic, AD_Statistic, GOF,
censored and at_boundary; IC and
KS_Statistic are kept.
Behaviour change. For censored data the Kolmogorov-Smirnov value changes: the supremum is now taken at both ends of every Kaplan-Meier step, where 0.3.0 looked only at the right-hand end and could understate it.
bd_analyze_csv() adds the two new statistics to its
goodness-of-fit table and a flag column to its
information-criteria table. compare_distributions() lists a
model without an AIC last and does not rank it.
bd_profile_ci() extends its default grid outward, in
the same logarithmic steps and starting each new point from its
neighbour, until the profile crosses the cut-off or the parameter box is
reached. Crossings are located by interpolation on the log scale and
each endpoint is flagged lower_closed /
upper_closed. New arguments extend and
max_extend. A grid supplied by the user is read exactly as
before. The interval is measured from the highest log-likelihood
reached; if the profile rises above the fitted maximum, a warning says
the fit stopped short of the optimum (new component
logLik_fit).
A standard error that is zero or comes from a degenerate
information matrix is reported as NA, never as
0, in summary(), bd_wald_ci(),
bd_identified_coef(), the AFT and cure summaries and
bd_analyze_csv().
Behaviour change.
bayes_betadanish() now uses, by default, the multi-chain
adaptive Metropolis sampler of the thesis
(method = "adaptive", four chains), written in base R. The
proposal scale adapts towards an acceptance rate of 0.234 and the
proposal covariance is learned during the burn-in; adaptation stops
before any draw is retained. Every fit reports R-hat, the effective
sample size and the acceptance rate of each parameter, and
converged is TRUE only when every R-hat is at
most 1.05 and every effective sample size at least 400. For the full
model the summary includes the identified composite
a_times_c. The MCMCpack sampler of earlier versions remains
available as method = "metrop1R". New arguments
method, n_chains, thin,
target_acc, rhat_max and ess_min;
the existing arguments keep their positions.
plot.bd_bayes() draws each chain in its own colour
in the trace panels.
bd_analyze_csv(bayes = TRUE) no longer requires
MCMCpack.
The Bayesian vignette is rewritten for the new sampler.
pbetadanish() stays finite where the Danish baseline
underflows (log G < -700). It uses the leading term of
the convergent expansion of the incomplete beta function there, so the
lower tail and its logarithm no longer collapse to 0 and
-Inf. Values elsewhere are unchanged. An undefined argument
(for example c and k both infinite) gives
NaN, as before, rather than an error.
qbetadanish() checks that lower.tail
and log.p are single TRUE/FALSE
values.
bd_kl_submodel() computes the Kullback-Leibler
divergence from a Beta-Danish distribution to its closest ED submodel by
numerical integration, and the sample size at which the four-parameter
model would be expected to win by a given AIC margin.The package now cites a single reference, the thesis (Ahmad,
2026), in DESCRIPTION, citation("BetaDanish"),
the help pages and the README.
The README lists the six vignettes that actually ship.
bd_ttt_plot() draws the scaled total time on test
transform, a distribution-free read on hazard shape before any model is
fitted. It returns the plotting coordinates and a suggested shape label,
so it can be compared against bd_hazard_shape() from the
fitted parameters.
bd_profile_plot() draws a profile log-likelihood
with the critical threshold and the resulting interval marked. Where the
curve never returns below the threshold the plot says there is no finite
upper bound, rather than implying one at the edge of the grid.
plot() method for bd_bayes objects,
giving trace and posterior density panels and restoring the graphical
parameters on exit.
The introduction and competing-risks vignettes cover the new functions.
fit_bd_competing() accepts covariates, either as a
one-sided formula with data or as a matrix. Each cause
carries its own coefficient vector acting on the time scale,
T_j = T_0j * exp(x'gamma_j), so a covariate may accelerate
one cause and retard another. The intercept is dropped, since each cause
already has its own scale parameter.
fit_bd_competing(submodel = TRUE) fits Exponentiated
Danish cause-specific kernels, and the fit now goes through the same
degeneracy guard as the univariate model.
cif_betadanish() gains an x argument
for the covariate values at which to evaluate the cumulative incidence,
defaulting to the reference subject. A covariate effect folds exactly
into the scale parameter, since scaling time by lambda maps
k to k/lambda.
Simulation studies:
bd_simulation_study() for the univariate model,
bd_simulation_cure() and
bd_simulation_competing(), alongside
simulate_bd_competing_data(). Each reports bias, RMSE, the
ratio of mean estimated standard error to empirical standard deviation,
and Wald coverage. Replicates with no admissible optimum are counted
rather than dropped.
fit_betadanish(grouped = TRUE) uses the grouped
likelihood for times recorded on a coarse grid. An event at
t contributes
log{F(t + delta/2) - F(t - delta/2)} rather than
log f(t). The increment is inferred from the spacing of the
times unless delta is supplied, and a warning now fires
when the times look grid-recorded but the point-density likelihood is
being used, which understates every standard error.
fit_betadanish(penalty = ) adds a ridge penalty on
the log-parameter scale for the weakly identified regime, shrinking
toward the unpenalised optimum unless penalty_center says
otherwise.
The reported logLik is the unpenalised
log-likelihood evaluated at the penalised estimate, so AIC, BIC and
likelihood ratio tests stay comparable with unpenalised fits. The
objective actually maximised is kept separately as
penalised_logLik.
bd_wald_ci() builds Wald intervals on the log scale
and exponentiates, so they cannot cross zero. The symmetric interval on
the natural scale does exactly that for a weakly identified tail
index.
bd_profile_ci() inverts the likelihood ratio test
for one parameter. When the profile stays inside the critical region to
the top of the grid it returns Inf and says so: an
unbounded interval is a result, and the honest report is a lower bound
rather than a point estimate.
bd_identified_coef() reports the four-parameter fit
as (ac, b, k). Near the lower tail a and
c enter almost entirely through their product, so the
composite is estimated precisely while the two separately are not. The
log-likelihood is unchanged: this is a reparametrisation, not a
different model.
Moments: bd_moments(),
bd_moment_summary(), bd_incomplete_moment()
and bd_conditional_moment(). Every integral is taken on the
finite Beta(a, b) scale under the substitution u = G(z), so
there is no series truncation. The existence condition
b > r is enforced: a moment that does not exist returns
Inf rather than a large finite number from a truncated
sum.
Ageing and inequality: bd_mrl() and
bd_rmrl() for mean residual life and mean inactivity time,
bd_mean_deviation(), bd_lorenz(),
bd_bonferroni() and bd_pwm().
Entropies: bd_entropy_shannon() now
uses the closed form, with method = "quadrature" retained
as an independent cross-check. bd_entropy_renyi() and
bd_entropy_tsallis() are new.
The quadrature bd_entropy_shannon() that shipped in
0.2.0 from R/entropy.R is superseded by the closed form and
that file is removed. The signature changes from
(a, b, c, k, subdivisions, rel.tol) to
(a, b, c, k, terms, method, rel.tol, subdivisions). Named
calls are unaffected; a positional call passing
subdivisions fifth would now set terms. The
old behaviour is available as method = "quadrature".
The closed form’s series terms decay like i^-(b+1), so
plain truncation is accurate for large b and not for small
b. Against high-precision integration the untruncated-tail
sum at M = 2000 is out by about 1e-8 at
b = 3 but by about 0.11 at
b = 0.5. An analytic Euler-Maclaurin tail is now added,
restoring agreement to roughly eight digits throughout.
Reliability: bd_stress_strength()
for R = P(X > Y) with X the strength,
evaluated on the Beta scale of the stress variable.
Shape: bd_hazard_shape() classifies
the hazard as increasing, decreasing, bathtub or upside-down bathtub
using Glaser’s eta(t) = -f'(t)/f(t), formed analytically,
and reports the mode.
Order statistics:
bd_order_stat_cdf() and
bd_order_stat_moments(), alongside the existing
bd_order_stat_pdf().
Tail behaviour: bd_tail_index()
records that the survival function is regularly varying with index
-b, that E(Z^r) is finite exactly when
b > r, that the moment generating function does not
exist, and that the distribution lies in the Frechet domain of
attraction.
Exponentiated Danish API: ded(),
ped(), qed(), red(),
sed() and hed() name the a = 1
submodel directly instead of requiring
dbetadanish(x, 1, b, c, k).
brain_cancer has been removed. The
dataset of 500 brain cancer patients that shipped in 0.1.0 and 0.2.0 is
no longer included, and the AFT example in the README now uses
melanoma. The entry under 0.1.0 below is left as written,
since it records what that release actually contained.
guinea_pig added (Bjerkedal 1960, n
= 72): survival times in days of guinea pigs injected with virulent
tubercle bacilli. This is the reference case for a well-identified
four-parameter fit, with (SE 1.20) sitting about 2.2 standard errors
clear of the ridge. On remission and
carbon_fibres the parent model sits on the flat direction
instead, so without this dataset no built-in example showed the full
model behaving well.
inst/extdata/ added with four
example CSVs – complete, censored, covariate and competing-risks – so
that examples and vignettes can demonstrate the file-driven
workflow.
bd_analyze_csv() reads a delimited
or Excel file, fits the requested model, assembles tidy result tables,
and optionally writes those tables and the diagnostic figures to a
directory. It covers the univariate, AFT, cure and competing-risks
paths. Nothing is written to disk unless output_dir is
supplied, and each fit is attempted independently so one failure is
recorded rather than fatal.
bd_csv_template() writes a
correctly shaped skeleton CSV in any of four layouts, so the expected
column structure can be seen rather than read about.
print, summary and plot
methods for the new bd_analysis class.
New vignette: “Analysing Your Own Data from a CSV File”.
read_survival_data() gains cause_col,
sep, dec, encoding,
na.strings, drop_na and quiet,
and will guess the time and status columns when they are not named.
Guessing deliberately excludes names such as censor, whose
polarity differs between conventions.
The bd_data_report attribute now also records which
columns were used and whether the times look recorded on a coarse
grid.
read_survival_data(covar_cols = "all") retains every
column that is not part of the response. Columns that are dropped are
now named in a message and recorded in the report, rather than
disappearing silently.
Degenerate optima are no longer reported as
fits. The Beta-Danish likelihood has ridges along which the
objective can be driven arbitrarily high with no finite maximiser: as
k and c grow with c/k fixed, the
family collapses to a Frechet-type limit.
optim_multistart() selected purely on the highest
log-likelihood behind nothing but a finiteness test, so a runaway beat
every honest optimum. A four-parameter fit on the
transplant data returned a log-likelihood of
+5.2e77 on 91 observations, and the resulting likelihood
ratio test reported a chi-squared statistic of 1.0e78.
fit_betadanish(), fit_bd_aft() and
fit_bd_cure() now search over a deterministic start grid
before any random starts, and accept only admissible optima: shape
parameters bounded by 500, scale within 1e-8 to
1e8, and a plausible per-observation log-likelihood. The
thresholds and the start grid are those already used by the thesis
analysis pipeline, so the two now agree.
The fitted object records starts_ok,
starts_rejected and loglik_spread. Warnings
fire when starts are discarded as degenerate and when the accepted
optima span more than two log-likelihood units, which is the signal that
n_starts is too low for the surface.
plot.bd_aft() and plot.bd_cure() now
produce the Cox-Snell residual plot. Both looked up the shape parameters
as coefficients["b"] and coefficients["c"],
but fit_bd_aft() and fit_bd_cure() store them
as log_b and log_c. The lookup returned
NA, the internal guard caught it, and the functions
returned without drawing anything. The lookup is fixed and the values
are exponentiated back to the natural scale.
dbetadanish() is now accurate in the far right tail.
The term was formed as log1p(-exp(log_G)), which pins at
about log(1e-16) once G rounds to one,
flooring the log-density near -36.8 regardless of its true value. It is
now formed as log(-expm1(-c * log1p(1/(k*t)))), which holds
full relative precision for every k*t > 0.
pbetadanish() computes the survival function through
the beta mirror identity , so a probability near one is never subtracted
from one.
qbetadanish() obtains directly from with the tail
flag reversed. The previous route computed
y^(1/c) / (k * (1 - y^(1/c))), which loses all significant
digits as u approaches one. Round-trip accuracy now holds
to p = 1 - 1e-13.
hbetadanish() no longer substitutes
-700 for an exhausted log-survival. A positive density with
zero survival gives Inf, an honest divergent hazard, and
the indeterminate case gives NaN.
The density, distribution, quantile, survival and hazard
functions now recycle their parameters element-wise against
x, q or p, following the usual
convention for R distribution functions.
fit_betadanish() records npar,
nobs, nevent, AIC and
BIC on the fitted object, and
report_betadanish() computes information criteria through
the logLik method. Both previously omitted AIC and BIC
silently.
summary.bd_aft() and summary.bd_cure()
report shape parameters on the natural scale with delta-method standard
errors and an exponentiated log-scale confidence interval, separately
from the regression coefficients. No Wald test is reported for a shape
parameter, since the implied null lies outside the parameter
space.
read_survival_data() selects columns by name rather
than by position, so a covariate named time or
status no longer collides with the response. Covariates
keep their original type. The returned data frame carries a
bd_data_report attribute.
extract_surv_data() validates times and the event
indicator up front and reports a clear error rather than failing inside
the optimiser.
sbetadanish(), an explicit survival
function.
fit_betadanish() gains
check_identifiability. It warns when the optimiser reports
a poor code, when the information matrix is singular, when is within two
standard errors of the non-identifiability ridge, and when the fitted
correlation between and exceeds 0.95 in absolute value. The
Identifiability section of ?fit_betadanish explains each
case.
?fit_betadanish documents the ridge and the
lower-tail confounding.
?fit_bd_competing documents the Tsiatis (1975)
non-identifiability result and the direction of bias under latent
dependence.
?BetaDanish records that in the upper tail and hence
that is finite if and only if .
The a = 1 submodel is called the Exponentiated
Danish (ED) throughout, in line with the underlying thesis.
“Complementary Exponentiated Danish (CED)” has been removed.
Carbon fibre breaking stress is given in GPa. It was previously
written “Gba”, and the typo was whitelisted in
inst/WORDLIST, which is why the spell check never caught
it.
The melanoma help page reported six columns and
documented seven.
Regression tests pin the tail behaviour of the distribution functions. The survival function is regularly varying with index , so the slope of against must approach . The previous implementation cannot pass this test, which is why the tail defects went unnoticed.
Smoke tests added for fit_bd_aft(),
fit_bd_cure(), plot.bd_aft() and
plot.bd_cure(), none of which were previously
covered.
Changelog correction. The 0.2.0 entry below has been rewritten to describe only what that release actually contained. As first published it also listed mean residual life, hazard-shape classification, stress-strength reliability, bootstrap confidence intervals for AFT and cure models, and a finite-sample simulation-study runner, none of which were implemented, and it reported four bug fixes that had not been applied. Those features are being added in the 0.3.0 development series and the fixes are recorded above.
bayes_betadanish()
provides random-walk Metropolis sampling for the Exponentiated Danish
submodel and the full four-parameter Beta-Danish model with vague Gamma
priors.fit_bd_competing() uses bound-constrained multi-start
L-BFGS-B optimisation. cif_compare() overlays fitted
cumulative incidence functions against the Aalen-Johansen estimator and
reports Gray’s test.bd_entropy_shannon(), by adaptive quadrature) and
order-statistic densities (bd_order_stat_pdf()).Three new vignettes were added:
Suggests packages guarded with
requireNamespace() at the call sites.