Semiparametric Bayesian Regression for Dependent Current Status Data

The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Introduction

The SemiParamBernsteinDepCS package provides a semiparametric Bayesian regression framework using Bernstein polynomial baseline models for analyzing dependent current status data. The methodology accommodates proportional hazards (PH) and proportional odds (PO) regression models paired with Archimedean copulas (Gumbel, Frank, and Clayton) to capture dependence between event time and observation/censoring time.

For full mathematical details, see Sharma and Balakrishnan (2026) doi:10.1080/02664763.2026.2701921.

Simulating Dependent Current Status Data

We generate a sample dataset under a Proportional Hazards (PH) model with a Gumbel copula:

library(SemiParamBernsteinDepCS)

set.seed(2026)
sim_data <- sim_bernstein_depcs(
  n = 100,
  model_type = "PH",
  copula = "gumbel",
  beta_Y = c(0.5, -0.3),
  beta_T = c(-0.2, 0.4)
)

head(sim_data)
#>     time status delta      x1      x2
#> 1 0.3647      1     0  0.4165  0.9730
#> 2 0.5764      1     0 -0.8638 -0.3951
#> 3 1.1980      1     1  0.1114 -0.1822
#> 4 0.2923      1     1 -0.0678 -0.8893
#> 5 0.4110      1     0 -0.5333  0.7967
#> 6 0.0836      1     0 -2.0129  0.4495

Model Fitting

We fit the semiparametric Bayesian regression model using fit_semiparam_bernstein_depcs():

fit_ph <- fit_semiparam_bernstein_depcs(
  formula_Y = delta ~ x1 + x2,
  formula_T = time ~ x1 + x2,
  data = sim_data,
  model_type = "PH",
  copula = "gumbel",
  order_m = 2,
  n_iter = 400,
  n_burn = 100,
  seed = 42
)

print(fit_ph)
#> Semiparametric Bayesian Regression Model for Dependent Current Status Data
#> =========================================================================
#> Model Type   : PH 
#> Copula Family: gumbel 
#> Bernstein Deg: Y=2, T=2 
#> Observations : 100 
#> MCMC Draws   : 300 (Burn-in: 100 , Thin: 1 )
#> Accept Rate  : 0.34 
#> DIC          : 251.743 
#> 
#> Posterior Means of Regression Parameters & Dependence:
#> beta_Y_x1 beta_Y_x2 beta_T_x1 beta_T_x2       S_Y       S_T     alpha       tau 
#>    0.1776   -0.0139   -0.1989    0.4990    5.1480    3.6114    3.7782    0.7332

Summary & Model Diagnostics

Extract posterior means, standard deviations, and 95% HPD credible intervals:

summary(fit_ph)
#> Posterior Summaries (95% HPD Credible Intervals):
#> =========================================================
#>              Mean     SD  Median HPD_lower_95 HPD_upper_95
#> beta_Y_x1  0.1776 0.1136  0.1838      -0.0398       0.3545
#> beta_Y_x2 -0.0139 0.1580 -0.0308      -0.3005       0.2788
#> beta_T_x1 -0.1989 0.1290 -0.2039      -0.3909       0.0125
#> beta_T_x2  0.4990 0.1401  0.4745       0.1915       0.7492
#> S_Y        5.1480 0.7118  5.1329       3.7276       6.6884
#> S_T        3.6114 0.6400  3.5443       2.2856       4.6717
#> alpha      3.7782 0.3405  3.7321       3.2362       4.4555
#> tau        0.7332 0.0242  0.7321       0.6910       0.7756
#> 
#> Model Selection Metrics:
#>   Log-Likelihood (at posterior mean): -119.041 
#>   Effective number of parameters pD : 6.831 
#>   DIC                               : 251.743

Prediction & Forecasting

Compute predicted marginal survival curves and HPD bounds for new covariate values:

pred_res <- predict(fit_ph, times = seq(0.1, 2.5, length.out = 20))
head(pred_res$predictions_Y)
#>        time      mean    median     lower     upper
#> 1 0.1000000 0.9596390 0.9608597 0.9409836 0.9728335
#> 2 0.2263158 0.8987311 0.8997490 0.8610672 0.9251443
#> 3 0.3526316 0.8290643 0.8298560 0.7766280 0.8647276
#> 4 0.4789474 0.7533395 0.7548100 0.6896208 0.8001646
#> 5 0.6052632 0.6743036 0.6753224 0.6091754 0.7394365
#> 6 0.7315789 0.5945808 0.5929988 0.5199709 0.6618729

Plotting

Visualize predicted marginal survival curves:

plot(fit_ph, type = "survival")

Real Data Example (Primary Biliary Cirrhosis)

The package includes the benchmark pbc_depcs dataset:

data(pbc_depcs)
head(pbc_depcs)
#>   time status delta age sex ascites bilirubin albumin
#> 1 3192      0     0  55   1       0      1.17    3.44
#> 2 1214      0     0  39   0       0      2.38    3.43
#> 3 2513      0     0  51   0       0      1.37    3.16
#> 4 2801      1     1  49   1       0      0.37    3.88
#> 5  594      0     1  43   0       1      1.31    3.16
#> 6 3235      1     0  25   0       0      6.78    3.20

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.