officially support the 'cox' family #230

NAMESPACE

@@ -323,6 +323,7 @@ export(cor_lagsar)
 export(cor_ma)
 export(cor_sar)
 export(cosy)
+export(cox)
 export(cratio)
 export(cs)
 export(cse)

NEWS.md

@@ -2,6 +2,8 @@
 
 ### New Features
 
+* Support the Cox proportional hazards model for 
+time-to-event data via family `cox`. (#230, #962)
 * Support method `loo_moment_match`, which can be used to
 update a `loo` object when Pareto k estimates are large.
 

R/families.R

@@ -82,9 +82,10 @@
 #'   \code{sratio} ('stopping ratio'), and \code{acat} ('adjacent category') 
 #'   leads to ordinal regression.} 
 #'   
-#'   \item{Families \code{Gamma}, \code{weibull}, \code{exponential}, 
-#'   \code{lognormal}, \code{frechet}, and \code{inverse.gaussian} can be 
-#'   used (among others) for survival regression.}
+#'   \item{Families \code{Gamma}, \code{weibull}, \code{exponential},
+#'   \code{lognormal}, \code{frechet}, \code{inverse.gaussian}, and \code{cox}
+#'   (Cox proportional hazards model) can be used (among others) for
+#'   time-to-event regression also known as survival regression.}
 #'   
 #'   \item{Families \code{weibull}, \code{frechet}, and \code{gen_extreme_value}
 #'   ('generalized extreme value') allow for modeling extremes.}
@@ -152,6 +153,9 @@
 #'   \item{Family \code{von_mises} supports \code{tan_half} and 
 #'   \code{identity}.}
 #'   
+#'   \item{Family \code{cox} supports \code{log}, \code{identity},
+#'   and \code{softplus} for the proportional hazards parameter.}
+#'   
 #'   \item{Family \code{wiener} supports \code{identity}, \code{log}, 
 #'   and \code{softplus} for the main parameter which represents the
 #'   drift rate.}
@@ -636,9 +640,8 @@ zero_inflated_asym_laplace <- function(link = "identity", link_sigma = "log",
               link_zi = link_zi)
 }
 
-# do not export yet!
-# @rdname brmsfamily
-# @export
+#' @rdname brmsfamily
+#' @export
 cox <- function(link = "log", bhaz = NULL) {
   slink <- substitute(link)
   .brmsfamily("cox", link = link, bhaz = bhaz)

R/family-lists.R

@@ -296,7 +296,7 @@
 
 .family_cox <- function() {
   list(
-    links = c("log", "identity"),
+    links = c("log", "identity", "softplus"),
     dpars = c("mu"), type = "real",
     ybounds = c(0, Inf), closed = c(TRUE, NA),
     ad = c("weights", "subset", "cens", "trunc"),

vignettes/brms_families.Rmd

@@ -94,11 +94,11 @@ $$
 with location parameter $\mu \in [0, 1]$ and positive shape parameter $\alpha$.
 -->
 
-## Survival models
+## Time-to-event models
 
-With survival models we mean all models that are defined on the positive reals
-only, that is $y \in \mathbb{R}^+$. The density of the **lognormal** family is
-given by
+With time-to-event models we mean all models that are defined on the positive
+reals only, that is $y \in \mathbb{R}^+$. The density of the **lognormal**
+family is given by
 $$
 f(y) = \frac{1}{\sqrt{2\pi}\sigma x} \exp\left(-\frac{1}{2}\left(\frac{\log(y) - \mu}{\sigma}\right)^2\right)
 $$ 
@@ -122,7 +122,15 @@ is set to $1$ for either the gamma or Weibull distribution. The density of the
 $$
 f(y) = \left(\frac{\alpha}{2 \pi y^3}\right)^{1/2} \exp \left(\frac{-\alpha (y - \mu)^2}{2 \mu^2 y} \right)
 $$
-where $\alpha$ is a positive shape parameter.
+where $\alpha$ is a positive shape parameter. The **cox** family implements Cox
+proportional hazards model which assumes a hazard function of the form $h(y) =
+h_0(y) \mu$ with baseline hazard $h_0(y)$ expressed via M-splines (which
+integrate to I-splines) in order to ensure monotonicity. The density of the cox
+model is then given by
+$$
+f(y) = h(y) S(y)
+$$
+where $S(y)$ is the survival function implied by $h(y)$.
 
 ## Extreme value models