Most users arrive with a spreadsheet rather than an R data frame. This vignette covers the file-driven path: from a CSV on disk to fitted models, tables and figures, without writing any modelling code.
Two columns are enough:
| column | meaning |
|---|---|
time |
survival or failure time, strictly positive |
status |
1 if the event was observed, 0 if right-censored |
Covariates go in additional columns. For competing risks, replace
status with cause, coded 0 for censored and 1,
2, … for the causes.
If you would rather start from a working example than a description, write a template and fill it in:
read_survival_data() will guess the time and status
columns when they carry recognisable names, and reports what it
chose:
f <- system.file("extdata", "censored_sample.csv", package = "BetaDanish")
dat <- read_survival_data(f)
head(dat)
#> time status
#> 1 11 1
#> 2 14 1
#> 3 23 1
#> 4 31 1
#> 5 32 1
#> 6 35 1Note that covariates are not retained by default. Ask for them by
name, or use covar_cols = "all":
dat_all <- read_survival_data(f, covar_cols = "all", quiet = TRUE)
names(dat_all)
#> [1] "time" "status" "group"The attached report records what happened, which is worth checking before you model anything:
str(attr(dat_all, "bd_data_report"))
#> List of 13
#> $ file : chr "C:\\Users\\hp\\AppData\\Local\\Temp\\RtmpSsnbRT\\Rinst50843b8c2feb\\BetaDanish\\extdata\\censored_sample.csv"
#> $ time_col : chr "time"
#> $ status_col : chr "status"
#> $ cause_col : NULL
#> $ available_columns: chr [1:3] "time" "status" "group"
#> $ covariates : chr "group"
#> $ dropped_columns : chr(0)
#> $ rows_read : int 91
#> $ rows_kept : int 91
#> $ rows_dropped : int 0
#> $ n_events : int 69
#> $ censoring_prop : num 0.242
#> $ grid_step : num 1One field deserves attention. grid_step is
non-NA when the times look recorded on a coarse grid –
whole days or whole months, say. Rounded times break the assumptions
behind a point-density likelihood, so treat standard errors from such
data as optimistic.
censorColumn-name guessing deliberately ignores names like
censor. In some conventions censor = 1 means
censored; in others it means observed. Guessing wrong
would silently invert every event in your dataset and still produce a
plausible-looking fit, so you must name that column yourself.
bd_analyze_csv() does the whole thing: read, fit,
tabulate, and optionally write results to a directory.
res <- bd_analyze_csv(f, analysis = "univariate", model = "ED",
compare = FALSE, n_starts = 5, seed = 1, quiet = TRUE)
res
#>
#> BetaDanish CSV analysis
#> -----------------------
#> Analysis: univariate
#> Observations: 91 (events: 69, censored: 24.2%)
#> Models fitted: ED
#>
#> model n npar logLik AIC BIC
#> ED 91 3 -476.8013 959.6026 967.1352
#>
#> Warnings: 1 (see $warnings)The tidy tables are in $tables:
res$tables$estimates
#> model parameter estimate std_error lower_95 upper_95
#> ED.1 ED b 0.5297814 0.08547109 0.3622611 0.6973016
#> ED.2 ED c 17.4808483 22.38401058 -26.3910063 61.3527029
#> ED.3 ED k 0.2597284 0.35870994 -0.4433301 0.9627870
res$tables$goodness_of_fit
#> model AIC BIC HQIC AICC KS
#> ED ED 959.6026 967.1352 962.6415 959.8784 0.0814148With model = "both", both the four-parameter Beta-Danish
and the three-parameter Exponentiated Danish submodel are fitted and
compared:
both <- bd_analyze_csv(f, model = "both", compare = FALSE,
n_starts = 5, seed = 2, quiet = TRUE)
#> Likelihood Ratio Test (a = 1 vs a != 1)
#>
#> Model LogLik Chisq Df Pr(>Chisq)
#> 1 Submodel (3-param) -476.8013 NA NA NA
#> 2 Full Model (4-param) -476.8225 0 1 1
both$tables$information_criteria
#> model n npar logLik AIC BIC
#> 1 BD 91 4 -476.8225 961.6450 971.6884
#> 2 ED 91 3 -476.8013 959.6026 967.1352
both$tables$likelihood_ratio_test
#> Model LogLik Chisq Df Pr(>Chisq)
#> 1 Submodel (3-param) -476.8013 NA NA NA
#> 2 Full Model (4-param) -476.8225 0 1 1Read that test alongside the identifiability diagnostics. On many
datasets the four-parameter model fits well but is only weakly
identified, and the likelihood ratio test will not reject the submodel.
See the Identifiability section of ?fit_betadanish, and the
guinea_pig dataset for a contrasting case where the parent
model is well identified.
Nothing is written to disk unless you name a directory. When you do, you get one CSV per table and a PNG per diagnostic figure:
out <- file.path(tempdir(), "bd_results")
saved <- bd_analyze_csv(f, model = "ED", compare = FALSE,
output_dir = out, n_starts = 5, seed = 3, quiet = TRUE)
basename(saved$files)
#> [1] "estimates.csv" "information_criteria.csv"
#> [3] "goodness_of_fit.csv" "data_report.csv"
#> [5] "survival.png" "hazard.png"
#> [7] "density.png" "pp.png"
#> [9] "qq.png"The same entry point handles the regression models. Covariates are taken from the file, either all of them or the ones you name:
g <- system.file("extdata", "covariate_sample.csv", package = "BetaDanish")
# Accelerated failure time
aft <- bd_analyze_csv(g, analysis = "aft",
covariates = c("age", "thickness"))
# Mixture cure model, cure fraction depending on ulceration
cure <- bd_analyze_csv(g, analysis = "cure", cure_formula = ~ ulcer)
# Competing risks
h <- system.file("extdata", "competing_sample.csv", package = "BetaDanish")
cr <- bd_analyze_csv(h, analysis = "competing", cause_col = "cause")