Semiparametric Survival Analysis Using Bernstein Polynomial
spbp(formula, ...)
| formula | a Surv object with time to event, status and explanatory terms. |
|---|---|
| ... | Arguments passed to `rstan::sampling` (e.g. iter, chains) or `rstan::optimizing`. |
An object of class 'spbp'.
Fits Bernstein Polynomial based Proportional regression to survival data.
library("spsurv") data("veteran") ## imports from survival package fit_mle <- spbp(Surv(time, status) ~ karno + factor(celltype), data = veteran, model = "po")#>summary(fit_mle)#> Bernstein Polynomial based Proportional Odds model #> Call: #> spbp.default(formula = Surv(time, status) ~ karno + factor(celltype), #> data = veteran, model = "po", approach = "mle") #> #> n= 137, number of events= 128 #> #> coef exp(coef) se(coef) z Pr(>|z|) #> karno -0.06127 0.94057 0.00877 -6.98 2.9e-12 *** #> factor(celltype)smallcell 1.28876 3.62828 0.43726 2.95 0.0032 ** #> factor(celltype)adeno 1.43812 4.21275 0.47370 3.04 0.0024 ** #> factor(celltype)large 0.10805 1.11411 0.46658 0.23 0.8169 #> --- #> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 #> #> Likelihood ratio test= 73.3 on 4 df, p=5e-15 #> Wald test = 65.4 on 4 df, p=2e-13fit_bayes <- spbp(Surv(time, status) ~ karno + factor(celltype), data = veteran, model = "po", approach = "bayes", cores = 1, iter = 300, chains = 1, priors = list(beta = c("normal(0,4)"), gamma = "lognormal(0,4)"))#> #> SAMPLING FOR MODEL 'spbp' NOW (CHAIN 1). #> Chain 1: #> Chain 1: Gradient evaluation took 5.7e-05 seconds #> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.57 seconds. #> Chain 1: Adjust your expectations accordingly! #> Chain 1: #> Chain 1: #> Chain 1: Iteration: 1 / 300 [ 0%] (Warmup) #> Chain 1: Iteration: 30 / 300 [ 10%] (Warmup) #> Chain 1: Iteration: 60 / 300 [ 20%] (Warmup) #> Chain 1: Iteration: 90 / 300 [ 30%] (Warmup) #> Chain 1: Iteration: 120 / 300 [ 40%] (Warmup) #> Chain 1: Iteration: 150 / 300 [ 50%] (Warmup) #> Chain 1: Iteration: 151 / 300 [ 50%] (Sampling) #> Chain 1: Iteration: 180 / 300 [ 60%] (Sampling) #> Chain 1: Iteration: 210 / 300 [ 70%] (Sampling) #> Chain 1: Iteration: 240 / 300 [ 80%] (Sampling) #> Chain 1: Iteration: 270 / 300 [ 90%] (Sampling) #> Chain 1: Iteration: 300 / 300 [100%] (Sampling) #> Chain 1: #> Chain 1: Elapsed Time: 0.652766 seconds (Warm-up) #> Chain 1: 0.29776 seconds (Sampling) #> Chain 1: 0.950526 seconds (Total) #> Chain 1:#> Warning: The largest R-hat is 1.06, indicating chains have not mixed. #> Running the chains for more iterations may help. See #> http://mc-stan.org/misc/warnings.html#r-hat#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. #> Running the chains for more iterations may help. See #> http://mc-stan.org/misc/warnings.html#bulk-ess#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable. #> Running the chains for more iterations may help. See #> http://mc-stan.org/misc/warnings.html#tail-ess#> Warning: Relative effective sample sizes ('r_eff' argument) not specified. #> For models fit with MCMC, the reported PSIS effective sample sizes and #> MCSE estimates will be over-optimistic.#> Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.#> Warning: #> 4 (2.9%) p_waic estimates greater than 0.4. We recommend trying loo instead.summary(fit_bayes)#> Bayesian Bernstein Polynomial based Proportional Odds model #> #> Call: #> spbp.default(formula = Surv(time, status) ~ karno + factor(celltype), #> data = veteran, approach = "bayes", model = "po", priors = list(beta = c("normal(0,4)"), #> gamma = "lognormal(0,4)"), cores = 1, iter = 300, chains = 1) #> #> n= 137, number of events= 128 #> #> mode mean se_mean sd 50% n_eff Rhat #> karno -0.0601 -0.0626 0.00123 0.0124 -0.0617 101 1.011 #> factor(celltype)smallcell 1.1377 1.2198 0.04020 0.4419 1.2048 121 0.996 #> factor(celltype)adeno 1.2891 1.3943 0.04578 0.4806 1.3565 110 1.005 #> factor(celltype)large 0.0460 0.0362 0.04366 0.4865 0.0323 124 1.007 #> lowerHPD upperHPD #> karno -0.0873 -0.0401 #> factor(celltype)smallcell 0.4751 2.0376 #> factor(celltype)adeno 0.5628 2.3613 #> factor(celltype)large -0.9890 0.7972 #> --- #> mean_exp median_exp sd_exp lowerHPD_exp upperHPD_exp #> karno 0.939 0.94 0.0117 0.916 0.961 #> factor(celltype)smallcell 3.739 3.34 1.7996 1.280 7.028 #> factor(celltype)adeno 4.524 3.88 2.2545 1.235 9.351 #> factor(celltype)large 1.172 1.03 0.6607 0.372 2.219 #> --- #> #> Computed from 150 by 137 log-likelihood matrix #> #> Estimate SE #> elpd_waic -719.8 20.4 #> p_waic 10.6 1.3 #> waic 1439.7 40.9 #> #> 4 (2.9%) p_waic estimates greater than 0.4. We recommend trying loo instead. #> #> Computed from 150 by 137 log-likelihood matrix #> #> Estimate SE #> elpd_loo -719.9 20.4 #> p_loo 10.7 1.3 #> looic 1439.8 40.9 #> ------ #> Monte Carlo SE of elpd_loo is 0.3. #> #> Pareto k diagnostic values: #> Count Pct. Min. n_eff #> (-Inf, 0.5] (good) 133 97.1% 88 #> (0.5, 0.7] (ok) 4 2.9% 39 #> (0.7, 1] (bad) 0 0.0% <NA> #> (1, Inf) (very bad) 0 0.0% <NA> #> #> All Pareto k estimates are ok (k < 0.7). #> See help('pareto-k-diagnostic') for details.