Feature/survival

DESCRIPTION

@@ -10,6 +10,7 @@ Authors@R: c(person("Jonah", "Gabry", email = "[email protected]", role = "au
              person(given = "Jacqueline Buros", family = "Novik",
                     role = "ctb", comment = "R/stan_jm.R"),
              person("AstraZeneca", role = "ctb", comment = "R/stan_jm.R"),
+             person("Eren", "Elci", role = "ctb", comment = "R/stan_surv.R"),
              person("Trustees of", "Columbia University", role = "cph"),
              person("Simon", "Wood", role = "cph", comment = "R/stan_gamm4.R"),
              person("R Core", "Deveopment Team", role = "cph",
@@ -45,6 +46,8 @@ Imports:
     rstan (>= 2.32.0),
     rstantools (>= 2.1.0),
     shinystan (>= 2.3.0),
+    splines,
+    splines2 (>= 0.2.7),
     stats,
     survival (>= 2.40.1),
     RcppParallel (>= 5.0.1),
@@ -61,6 +64,7 @@ Suggests:
     mgcv (>= 1.8-13),
     rmarkdown,
     roxygen2,
+    simsurv (>= 0.2.2),
     StanHeaders (>= 2.21.0),
     testthat (>= 1.0.2),
     gamm4,

NAMESPACE

@@ -23,6 +23,8 @@ S3method(fixef,stanmvreg)
 S3method(fixef,stanreg)
 S3method(formula,stanmvreg)
 S3method(formula,stanreg)
+S3method(get_surv,stanjm)
+S3method(get_surv,stansurv)
 S3method(get_x,default)
 S3method(get_x,gamm4)
 S3method(get_x,lmerMod)
@@ -56,12 +58,16 @@ S3method(nsamples,stanreg)
 S3method(pairs,stanreg)
 S3method(plot,predict.stanjm)
 S3method(plot,stanreg)
+S3method(plot,stansurv)
 S3method(plot,survfit.stanjm)
+S3method(plot,survfit.stansurv)
 S3method(posterior_epred,stanreg)
 S3method(posterior_interval,stanreg)
 S3method(posterior_linpred,stanreg)
 S3method(posterior_predict,stanmvreg)
 S3method(posterior_predict,stanreg)
+S3method(posterior_survfit,stanjm)
+S3method(posterior_survfit,stansurv)
 S3method(posterior_vs_prior,stanreg)
 S3method(pp_check,stanreg)
 S3method(predict,stanreg)
@@ -79,6 +85,7 @@ S3method(print,stanreg_list)
 S3method(print,summary.stanmvreg)
 S3method(print,summary.stanreg)
 S3method(print,survfit.stanjm)
+S3method(print,survfit.stansurv)
 S3method(prior_summary,stanreg)
 S3method(ranef,stanmvreg)
 S3method(ranef,stanreg)
@@ -115,6 +122,7 @@ export(default_prior_intercept)
 export(dirichlet)
 export(exponential)
 export(fixef)
+export(get_surv)
 export(get_x)
 export(get_y)
 export(get_z)
@@ -182,10 +190,12 @@ export(stan_mvmer)
 export(stan_nlmer)
 export(stan_polr)
 export(stan_polr.fit)
+export(stan_surv)
 export(stanjm_list)
 export(stanmvreg_list)
 export(stanreg_list)
 export(student_t)
+export(tve)
 export(waic)
 if(getRversion()>='3.3.0') importFrom(stats, sigma) else
   importFrom(lme4,sigma)
@@ -195,6 +205,7 @@ import(bayesplot)
 import(methods)
 import(rstantools)
 import(shinystan)
+import(splines2)
 import(stats)
 importFrom(Matrix,Matrix)
 importFrom(Matrix,t)
@@ -272,6 +283,7 @@ importFrom(rstan,stanc)
 importFrom(rstan,vb)
 importFrom(rstantools,loo_R2)
 importFrom(rstantools,nsamples)
+importFrom(splines,bs)
 importFrom(stats,cov2cor)
 importFrom(stats,getInitial)
 importFrom(survival,Surv)

R/data_block.R

@@ -78,7 +78,8 @@ handle_glm_prior <- function(prior, nvars, default_scale, link,
                 prior_scale = as.array(rep(1, nvars)),
                 prior_df = as.array(rep(1, nvars)), prior_dist_name = NA,
                 global_prior_scale = 0, global_prior_df = 0,
-                slab_df = 0, slab_scale = 0,
+                slab_df = 0, slab_scale = 0, 
+                prior_concentration = as.array(rep(1, nvars)),
                 prior_autoscale = FALSE))
 
   if (!is.list(prior)) 
@@ -94,6 +95,7 @@ handle_glm_prior <- function(prior, nvars, default_scale, link,
   global_prior_df <- 0
   slab_df <- 0
   slab_scale <- 0
+  prior_concentration <- 1
   if (!prior_dist_name %in% unlist(ok_dists)) {
     stop("The prior distribution should be one of ",
          paste(names(ok_dists), collapse = ", "))
@@ -114,14 +116,24 @@ handle_glm_prior <- function(prior, nvars, default_scale, link,
     slab_scale <- prior$slab_scale
   } else if (prior_dist_name %in% "exponential") {
     prior_dist <- 3L # only used for scale parameters so 3 not a conflict with 3 for hs
+  } else if (prior_dist_name %in% "dirichlet") {
+    prior_dist <- 4L # only used by stan_surv for baseline hazard coefficients
+    prior_concentration <- prior$concentration
+    if (is.null(prior_concentration)) {
+      prior_concentration <- rep(1, nvars)
+    } else {
+      prior_concentration <- maybe_broadcast(prior_concentration, nvars)
+    }
   }
 
   prior_df <- maybe_broadcast(prior_df, nvars)
   prior_df <- as.array(pmin(.Machine$double.xmax, prior_df))
   prior_mean <- maybe_broadcast(prior_mean, nvars)
   prior_mean <- as.array(prior_mean)
   prior_scale <- maybe_broadcast(prior_scale, nvars)
-
+  prior_concentration <- maybe_broadcast(prior_concentration, nvars) 
+  prior_concentration <- as.array(prior_concentration)
+
   nlist(prior_dist, 
         prior_mean, 
         prior_scale, 
@@ -131,5 +143,6 @@ handle_glm_prior <- function(prior, nvars, default_scale, link,
         global_prior_df,
         slab_df,
         slab_scale,
+        prior_concentration,
         prior_autoscale = isTRUE(prior$autoscale))
 }

R/doc-datasets.R

@@ -20,7 +20,7 @@
 #' Small datasets for use in \pkg{rstanarm} examples and vignettes.
 #'
 #' @name rstanarm-datasets
-#' @aliases kidiq roaches wells bball1970 bball2006 mortality tumors radon pbcLong pbcSurv
+#' @aliases bball1970 bball2006 bcancer frail kidiq mice mortality pbcLong pbcSurv tumors radon roaches wells
 #' @format 
 #' \describe{
 #' \item{\code{bball1970}}{
@@ -50,6 +50,42 @@
 #' \item \code{K} Number of at-bats
 #' }
 #' }
+#' \item{\code{bcancer}}{
+#' The German Breast Cancer Study Group dataset, containing time to death or
+#' recurrence for 686 patients with primary node positive breast cancer 
+#' recruited between 1984-1989. 
+#' 
+#' Source: Royston and Parmar (2002)
+#' 
+#' 686 obs. of 4 variables
+#' \itemize{
+#' \item \code{recdays} Time to death or censoring (in days)
+#' \item \code{recyrs} Time to death or censoring (in years)
+#' \item \code{status} Event indicator (0 = right censored, 1 = event)
+#' \item \code{group} Prognostic group, based on a regression model developed
+#' by Sauerbrei and Royston (1999) (\code{Good}, \code{Medium}, \code{Poor})
+#' }
+#' }
+#' \item{\code{frail}}{
+#' A simulated dataset of event times (i.e. survival data) for 200 patients
+#' clustered within 20 hospital sites (10 patients per hospital site). 
+#' The event times are simulated from a parametric proportional hazards model 
+#' under the following assumptions: (i) a constant (i.e. exponential) baseline 
+#' hazard rate of 0.1; (ii) a fixed treatment effect with log hazard ratio of 
+#' 0.3; and (iii) a site-specific random intercept (specified on the log 
+#' hazard scale) drawn from a \eqn{N(0,1)} distribution.
+#' 
+#' 200 obs. of 6 variables
+#' \itemize{
+#' \item \code{id} ID unique to each patient
+#' \item \code{site} ID unique to each hospital site (i.e. cluster)
+#' \item \code{trt} Treatment indicator (0 = untreated, 1 = treated)
+#' \item \code{b} Cluster-specific random intercept used to simulate the 
+#' event times
+#' \item \code{eventtime} Event or censoring time
+#' \item \code{status} Event indicator (0 = right censored, 1 = event)
+#' }
+#' }
 #' \item{\code{kidiq}}{
 #' Data from a survey of adult American women and their children 
 #' (a subsample from the National Longitudinal Survey of Youth).
@@ -64,6 +100,24 @@
 #' \item \code{mom_age} Mother's age
 #' }
 #' }
+#' \item{\code{mice}}{
+#' Lung tumor development in 144 RFM mice allocated to either a conventional
+#' environment or germ-free environment. Mice were sacrificed and examined 
+#' for presence of a lung tumor. The outcome variables in the dataset 
+#' (\code{l} and \code{u}) denote a left-censored or right-censored time 
+#' interval within which the development of the first lung tumor must have 
+#' occurred. 
+#' 
+#' Source: Hoel and Walburg (1972)
+#'
+#' 144 obs. of 3 variables
+#' \itemize{
+#' \item \code{l} Lower limit of the interval.
+#' \item \code{u} Upper limit of the interval.
+#' \item \code{grp} Experimental group (\code{ce} = conventional environment,
+#' \code{ge} = germ-free environment).
+#' }
+#' }
 #' \item{\code{mortality}}{
 #' Surgical mortality rates in 12 hospitals performing cardiac surgery
 #' in babies.
@@ -86,21 +140,20 @@
 #' 
 #' 304 obs. of 8 variables (\code{pbcLong}) and 40 obs. of 7 variables (\code{pbcSurv})
 #' \itemize{
-#' \item \code{age} in years
-#' \item \code{albumin} serum albumin (g/dl)
-#' \item \code{logBili} logarithm of serum bilirubin
-#' \item \code{death} indicator of death at endpoint
-#' \item \code{futimeYears} time (in years) between baseline and  
+#' \item \code{age} Age (in years)
+#' \item \code{albumin} Serum albumin (g/dl)
+#' \item \code{logBili} Logarithm of serum bilirubin
+#' \item \code{death} Indicator of death at endpoint
+#' \item \code{futimeYears} Time (in years) between baseline and  
 #'     the earliest of death, transplantion or censoring
-#' \item \code{id} numeric ID unique to each individual
-#' \item \code{platelet} platelet count
-#' \item \code{sex} gender (m = male, f = female)
-#' \item \code{status} status at endpoint (0 = censored, 
-#'     1 = transplant, 2 = dead)
-#' \item \code{trt} binary treatment code (0 = placebo, 1 = 
-#'     D-penicillamine)
-#' \item \code{year} time (in years) of the longitudinal measurements,
-#'     taken as time since baseline)
+#' \item \code{id} Numeric ID unique to each individual
+#' \item \code{platelet} Platelet count
+#' \item \code{sex} Gender (m = male, f = female)
+#' \item \code{status} Status at endpoint (0 = censored, 1 = transplant, 
+#'     2 = dead)
+#' \item \code{trt} Binary treatment code (0 = placebo, 1 = D-penicillamine)
+#' \item \code{year} Time (in years) of the longitudinal measurements,
+#'     taken as time since baseline
 #' }
 #' }
 #' 
@@ -181,7 +234,20 @@
 #' @template reference-gelman-hill
 #' 
 #' @references
-#' Spiegelhalter, D., Thomas, A., Best, N., & Gilks, W. (1996) BUGS 0.5 
+#' Hoel, D. and Walburg, H. (1972) Statistical analysis of survival experiments. 
+#' \emph{The Annals of Statistics} \strong{18}:1259--1294.
+#' 
+#' Royston, P. and Parmar, M. (2002) Flexible parametric proportional-hazards 
+#' and proportional-odds models for censored survival data, with application 
+#' to prognostic modelling and estimation of treatment effects. 
+#' \emph{Statistics in Medicine} \strong{21}(1):2175--2197.
+#' 
+#' Sauerbrei, W. and Royston, P. (1999) Building multivariable prognostic and 
+#' diagnostic models: transformation of the predictors using fractional 
+#' polynomials. \emph{Journal of the Royal Statistical Society, Series A}
+#' \strong{162}:71--94.
+#' 
+#' Spiegelhalter, D., Thomas, A., Best, N., and Gilks, W. (1996) BUGS 0.5 
 #' Examples. MRC Biostatistics Unit, Institute of Public health, Cambridge, UK.
 #' 
 #' Tarone, R. E. (1982) The use of historical control information in testing for

R/doc-modeling-functions.R

@@ -1,7 +1,7 @@
 #' Modeling functions available in \pkg{rstanarm}
-#' 
+#'
 #' @name available-models
-#' 
+#'
 #' @section Modeling functions:
 #' The model estimating functions are described in greater detail in their
 #' individual help pages and vignettes. Here we provide a very brief
@@ -33,9 +33,16 @@
 #'   appropriate estimates of uncertainty for models that consist of a mix of
 #'   common and group-specific parameters.
 #'  }
+#'  \item{\code{\link{stan_mvmer}}}{
+#'    A multivariate form of \code{\link{stan_glmer}}, whereby the user can
+#'    specify one or more submodels each consisting of a GLM with group-specific
+#'    terms. If more than one submodel is specified (i.e. there is more than one
+#'    outcome variable) then a dependence is induced by assuming that the
+#'    group-specific terms for each grouping factor are correlated across submodels.
+#'  }
 #'  \item{\code{\link{stan_nlmer}}}{
-#'   Similar to \code{\link[lme4]{nlmer}} in the \pkg{lme4} package for 
-#'   nonlinear "mixed-effects" models, but the group-specific coefficients 
+#'   Similar to \code{\link[lme4]{nlmer}} in the \pkg{lme4} package for
+#'   nonlinear "mixed-effects" models, but the group-specific coefficients
 #'   have flexible priors on their unknown covariance matrices.
 #'  }
 #'  \item{\code{\link{stan_gamm4}}}{
@@ -68,19 +75,22 @@
 #'    to \code{\link[survival]{clogit}} that allow \code{stan_clogit} to accept
 #'    group-specific terms as in \code{\link{stan_glmer}}.
 #'  }
-#'  \item{\code{\link{stan_mvmer}}}{
-#'    A multivariate form of \code{\link{stan_glmer}}, whereby the user can
-#'    specify one or more submodels each consisting of a GLM with group-specific 
-#'    terms. If more than one submodel is specified (i.e. there is more than one
-#'    outcome variable) then a dependence is induced by assuming that the
-#'    group-specific terms for each grouping factor are correlated across submodels. 
+#'  \item{\code{\link{stan_surv}}}{
+#'    Fits models to survival (i.e. time-to-event) data using either a hazard 
+#'    scale or accelerated failure time (AFT) formulation. The user can choose 
+#'    between either a flexible spline-based approximation for the baseline 
+#'    hazard or a standard parametric distributional assumption for the baseline 
+#'    hazard. Covariate effects can then be accommodated under proportional or 
+#'    non-proportional hazards assumptions (for models on the hazard scale) or 
+#'    using time-fixed or time-varying acceleration factors (for models on the 
+#'    AFT scale).
 #'  }
 #'  \item{\code{\link{stan_jm}}}{
-#'    Estimates shared parameter joint models for longitudinal and time-to-event 
-#'    (i.e. survival) data. The joint model can be univariate (i.e. one longitudinal 
-#'    outcome) or multivariate (i.e. more than one longitudinal outcome). A variety 
-#'    of parameterisations are available for linking the longitudinal and event 
-#'    processes (i.e. a variety of association structures).      
+#'    Fits shared parameter joint models for longitudinal and survival (i.e.
+#'    time-to-event) data. The joint model can be univariate (i.e. one longitudinal
+#'    outcome) or multivariate (i.e. more than one longitudinal outcome). A variety
+#'    of parameterisations are available for linking the longitudinal and event
+#'    processes (i.e. a variety of association structures).
 #'  }
 #' }
 #'