Mixture cure models (Weibull and log-normal)

.gitignore

@@ -10,3 +10,4 @@ tests/local/models_0.10.0.Rda
 tests/local/models_1.2.0.Rda
 tests/local/Rplots.pdf
 /Meta/
+*ipynb*

DESCRIPTION

@@ -16,7 +16,8 @@ Authors@R:
       person("Mattan S.", "Ben-Shachar", role = c("ctb")),
       person("Hayden", "Rabel", role = c("ctb")),
       person("Simon C.", "Mills", role = c("ctb")),
-      person("Stephen", "Wild", role = c("ctb")))
+      person("Stephen", "Wild", role = c("ctb")),
+      person("Yanan", "Long",  role = c("ctb")))
 Depends:
     R (>= 3.6.0),
     Rcpp (>= 0.12.0),

NAMESPACE

@@ -387,6 +387,8 @@ export(dhurdle_poisson)
 export(dinv_gaussian)
 export(dirichlet)
 export(dlogistic_normal)
+export(dmixcure_lognormal)
+export(dmixcure_weibull)
 export(dmulti_normal)
 export(dmulti_student_t)
 export(do_call)
@@ -466,6 +468,8 @@ export(marginal_smooths)
 export(mcmc_plot)
 export(me)
 export(mi)
+export(mixcure_lognormal)
+export(mixcure_weibull)
 export(mixture)
 export(mm)
 export(mmc)
@@ -498,6 +502,8 @@ export(phurdle_lognormal)
 export(phurdle_negbinomial)
 export(phurdle_poisson)
 export(pinv_gaussian)
+export(pmixcure_lognormal)
+export(pmixcure_weibull)
 export(post_prob)
 export(posterior_average)
 export(posterior_epred)

R/distributions.R

@@ -2049,6 +2049,125 @@ phurdle_lognormal <- function(q, mu, sigma, hu, lower.tail = TRUE,
   out
 }
 
+#' Mixture cure (mixcure) Distributions
+#'
+#' Density and distribution functions for mixcure distributions.
+#'
+#' @name Mixcure
+#'
+#' @inheritParams StudentT
+#' @param inc incidence probability
+#' @param mu location parameter
+#' @param shape shape parameter
+#' @param sigma,scale scale parameter
+#'
+#' @details
+#' The mixcure distribution is a subclass of survival/time-to-event models
+#' where a subpopulation of the subjects are posited to be ``cured''
+#' i.e. event-free. Mathematically, 
+NULL
+
+#' @rdname Mixcure
+#' @export
+dmixcure_lognormal <- function(x, mu, sigma, inc, log = FALSE) {
+  pars <- list(meanlog = mu, sdlog = sigma)
+  .dmixcure(x, "lnorm", inc, pars, log)
+}
+
+#' @rdname Mixcure
+#' @export
+pmixcure_lognormal <- function(q, mu, sigma, inc, lower.tail = TRUE, log.p = FALSE) {
+  pars <- list(meanlog = mu, sdlog = sigma)
+  .pmixcure(q, "lnorm", inc, pars, lower.tail, log.p)
+}
+
+#' @rdname Mixcure
+#' @export
+dmixcure_weibull <- function(x, shape, scale, inc, log = FALSE) {
+  pars <- list(shape = shape, scale = scale)
+  .dmixcure(x, "weibull", inc, pars, log)
+}
+
+#' @rdname Mixcure
+#' @export
+pmixcure_weibull <- function(q, shape, scale, inc, lower.tail = TRUE, log.p = FALSE) {
+  pars <- list(shape = shape, scale = scale)
+  .pmixcure(q, "weibull", inc, pars, lower.tail, log.p)
+}
+
+#' @rdname Mixcure
+#' @export
+dmixcure_loglogistic <- function(x, shape, scale, inc, log = FALSE) {
+  pars <- list(shape = shape, scale = scale)
+  .dmixcure(x, "llogis", inc, pars, log)
+}
+
+#' @rdname Mixcure
+#' @export
+pmixcure_loglogistic <- function(q, shape, scale, inc, lower.tail = TRUE, log.p = FALSE) {
+  pars <- list(shape = shape, scale = scale)
+  .pmixcure(q, "llogis", inc, pars, lower.tail, log.p)
+}
+
+# density of a mixcure distribution
+# @param dist name of the distribution
+# @param inc bernoulli incidence parameter
+# @param pars list of parameters passed to pdf
+.dmixcure <- function(x, dist, inc, pars, log) {
+  stopifnot(is.list(pars))
+  dist <- as_one_character(dist)
+  log <- as_one_logical(log)
+  args <- expand(dots = c(nlist(x, inc), pars))
+  x <- args$x
+  inc <- args$inc
+  pars <- args[names(pars)]
+  pdf <- paste0("d", dist)
+  # incidence part (not censored): pi(z) * f(t | x)
+  out <- log(inc) + do_call(pdf, c(list(x), pars, log = TRUE))
+  if (!log) {
+    out <- exp(out)
+  }
+  out
+}
+
+# CDF of a mixcure distribution
+# @param dist name of the distribution
+# @param inc bernoulli incidence parameter
+# @param pars list of parameters passed to pdf
+# @param lb lower bound of the conditional distribution
+# @param ub upper bound of the conditional distribution
+.pmixcure <- function(q, dist, inc, pars, lower.tail, log.p, lb = 0, ub = Inf) {
+  stopifnot(is.list(pars))
+  dist <- as_one_character(dist)
+  lower.tail <- as_one_logical(lower.tail)
+  log.p <- as_one_logical(log.p)
+  args <- expand(dots = c(nlist(q, inc), pars))
+  q <- args$q
+  inc <- args$inc
+  pars <- args[names(pars)]
+  cdf <- paste0("p", dist)
+  # compute log CCDF values
+  # latency part (right-censored): [1 - pi(z)] + pi(z) * S(t | x)
+  out <- matrixStats::logSumExp(c(
+    1, -inc,
+    inc + do_call(cdf, c(list(q), pars, lower.tail = FALSE, log.p = TRUE))
+  ))
+  # # take the limits of the distribution into account
+  out <- ifelse(q < lb, 0, out)
+  out <- ifelse(q > ub, -Inf, out)
+  if (lower.tail) {
+    out <- 1 - exp(out)
+    if (log.p) {
+      out <- log(out)
+    }
+  } else {
+    if (!log.p) {
+      out <- exp(out)
+    }
+  }
+  out
+}
+
 # density of the categorical distribution with the softmax transform
 # @param x positive integers not greater than ncat
 # @param eta the linear predictor (of length or ncol ncat)

R/families.R

@@ -35,6 +35,7 @@
 #' @param link_beta Link of auxiliary parameter \code{beta} if being predicted.
 #' @param link_zi Link of auxiliary parameter \code{zi} if being predicted.
 #' @param link_hu Link of auxiliary parameter \code{hu} if being predicted.
+#' @param link_inc Link of auxiliary parameter \code{inc} if being predicted.
 #' @param link_zoi Link of auxiliary parameter \code{zoi} if being predicted.
 #' @param link_coi Link of auxiliary parameter \code{coi} if being predicted.
 #' @param link_disc Link of auxiliary parameter \code{disc} if being predicted.
@@ -192,7 +193,8 @@ brmsfamily <- function(family, link = NULL, link_sigma = "log",
                        link_shape = "log", link_nu = "logm1",
                        link_phi = "log", link_kappa = "log",
                        link_beta = "log", link_zi = "logit",
-                       link_hu = "logit", link_zoi = "logit",
+                       link_hu = "logit", link_inc = "logit",
+                       link_zoi = "logit",
                        link_coi = "logit", link_disc = "log",
                        link_bs = "log", link_ndt = "log",
                        link_bias = "logit", link_xi = "log1p",
@@ -206,7 +208,7 @@ brmsfamily <- function(family, link = NULL, link_sigma = "log",
     link_sigma = link_sigma, link_shape = link_shape,
     link_nu = link_nu, link_phi = link_phi,
     link_kappa = link_kappa, link_beta = link_beta,
-    link_zi = link_zi, link_hu = link_hu,
+    link_zi = link_zi, link_hu = link_hu, link_inc = link_inc,
     link_zoi = link_zoi, link_coi = link_coi,
     link_disc = link_disc, link_bs = link_bs,
     link_ndt = link_ndt, link_bias = link_bias,
@@ -603,6 +605,13 @@ gen_extreme_value <- function(link = "identity", link_sigma = "log",
               link_sigma = link_sigma, link_xi = link_xi)
 }
 
+#' @rdname brmsfamily
+#' @export
+loglogistic <- function(link = "log", link_shape = "log") {
+  slink <- substitute(link)
+  .brmsfamily("loglogistic", link = link, link_shape = link_shape)
+}
+
 #' @rdname brmsfamily
 #' @export
 exgaussian <- function(link = "identity", link_sigma = "log",
@@ -736,6 +745,33 @@ hurdle_cumulative <- function(link = "logit", link_hu = "logit",
                threshold = threshold)
 }
 
+#' @rdname brmsfamily
+#' @export
+mixcure_lognormal <- function(link = "identity", link_sigma = "log",
+                              link_inc = "logit") {
+  slink <- substitute(link)
+  .brmsfamily("mixcure_lognormal", link = link, slink = slink,
+              link_sigma = link_sigma, link_inc = link_inc)
+}
+
+#' @rdname brmsfamily
+#' @export
+mixcure_weibull <- function(link = "log", link_shape = "log",
+                            link_inc = "logit") {
+  slink <- substitute(link)
+  .brmsfamily("mixcure_weibull", link = link, slink = slink,
+              link_shape = link_shape, link_inc = link_inc)
+}
+
+#' @rdname brmsfamily
+#' @export
+mixcure_loglogistic <- function(link = "identity", link_shape = "log",
+                              link_inc = "logit") {
+  slink <- substitute(link)
+  .brmsfamily("mixcure_loglogistic", link = link, slink = slink,
+              link_shape = link_shape, link_inc = link_inc)
+}
+
 #' @rdname brmsfamily
 #' @export
 zero_inflated_beta <- function(link = "logit", link_phi = "log",
@@ -1334,6 +1370,7 @@ links_dpars <- function(dpar) {
     beta = c("log", "identity", "softplus", "squareplus"),
     zi = c("logit", "identity"),
     hu = c("logit", "identity"),
+    inc = c("logit", "cloglog", "probit", "probit_approx"),
     zoi = c("logit", "identity"),
     coi = c("logit", "identity"),
     disc = c("log", "identity", "softplus", "squareplus"),
@@ -1872,10 +1909,10 @@ family_bounds.brmsterms <- function(x, ...) {
   pos_families <- c(
     "poisson", "negbinomial", "negbinomial2", "geometric",
     "gamma", "weibull", "exponential", "lognormal",
-    "frechet", "inverse.gaussian",
+    "frechet", "inverse.gaussian", "loglogistic",
     "hurdle_poisson", "hurdle_negbinomial", "hurdle_gamma",
-    "hurdle_lognormal", "zero_inflated_poisson",
-    "zero_inflated_negbinomial"
+    "hurdle_lognormal", "mixcure_lognormal", "mixcure_weibull",
+    "mixcure_loglogistic", "zero_inflated_poisson", "zero_inflated_negbinomial"
   )
   beta_families <- c("beta", "zero_inflated_beta", "zero_one_inflated_beta")
   ordinal_families <- c("cumulative", "cratio", "sratio", "acat")

R/family-lists.R

@@ -330,6 +330,17 @@
   )
 }
 
+.family_loglogistic <- function() {
+  list(
+    links = c("log", "identity"),
+    dpars = c("mu", "shape"), type = "real",
+    ybounds = c(0, Inf), closed = c(FALSE, NA),
+    ad = c("weights", "subset", "cens", "trunc", "mi", "index"),
+    include = "fun_loglogistic.stan",
+    specials = "logscale"
+  )
+}
+
 .family_von_mises <- function() {
   list(
     links = c("tan_half", "identity"),
@@ -498,6 +509,30 @@
   )
 }
 
+.family_mixcure_lognormal <- function() {
+  list(
+    links = c("identity", "inverse"),
+    dpars = c("mu", "sigma", "inc"), type = "real",
+    ybounds = c(0, Inf), closed = c(TRUE, NA),
+    ad = c("weights", "subset", "cens", "trunc", "index"),
+    include = "fun_mixcure_lognormal.stan",
+    specials = c("logscale", "sbi_inc_logit", "sbi_inc_logit_cdf"),
+    normalized = ""
+  )
+}
+
+.family_mixcure_weibull <- function() {
+  list(
+    links = c("log", "identity", "inverse", "softplus", "squareplus"),
+    dpars = c("mu", "shape", "inc"), type = "real",
+    ybounds = c(0, Inf), closed = c(TRUE, NA),
+    ad = c("weights", "subset", "cens", "trunc", "index"),
+    include = "fun_mixcure_weibull.stan",
+    specials = c("logscale", "sbi_inc_logit", "sbi_inc_logit_cdf"),
+    normalized = ""
+  )
+}
+
 .family_zero_inflated_poisson <- function() {
   list(
     links = c("log", "identity", "sqrt", "softplus", "squareplus"),

R/log_lik.R

@@ -586,6 +586,15 @@ log_lik_gen_extreme_value <- function(i, prep) {
   log_lik_weight(out, i = i, prep = prep)
 }
 
+log_lik_loglogistic <- function(i, prep) {
+  shape <- get_dpar(prep, "shape", i = i)
+  mu <- get_dpar(prep, "mu", i)
+  args <- list(shape = shape, scale = mu)
+  out <- log_lik_censor(dist = "llogis", args = args, i = i, prep = prep)
+  out <- log_lik_truncate(out, cdf = pllogis, args = args, i = i, prep = prep)
+  log_lik_weight(out, i = i, prep = prep)
+}
+
 log_lik_inverse.gaussian <- function(i, prep) {
   args <- list(mu = get_dpar(prep, "mu", i),
                shape = get_dpar(prep, "shape", i = i))
@@ -742,6 +751,44 @@ log_lik_hurdle_cumulative <- function(i, prep) {
   log_lik_weight(out, i = i, prep = prep)
 }
 
+log_lik_mixcure_lognormal <- function(i, prep) {
+  mu <- get_dpar(prep, "mu", i)
+  sigma <- get_dpar(prep, "sigma", i = i)
+  inc <- get_dpar(prep, "inc", i)
+  args <- nlist(mu = mu, sigma = sigma, inc = inc)
+  out <- log_lik_censor("mixcure_lognormal", args, i, prep)
+  out <- log_lik_truncate(out, pmixcure_lognormal, args, i, prep)
+  log_lik_weight(out, i = i, prep = prep)
+}
+
+log_lik_mixcure_weibull <- function(i, prep) {
+  shape <- get_dpar(prep, "shape", i = i)
+  scale <- get_dpar(prep, "mu", i = i) / gamma(1 + 1 / shape)
+  inc <- get_dpar(prep, "inc", i)
+  args <- list(shape = shape, scale = scale, inc = inc)
+  out <- log_lik_censor(
+    dist = "mixcure_weibull", args = args, i = i, prep = prep
+  )
+  out <- log_lik_truncate(
+    out, cdf = pmixcure_weibull, args = args, i = i, prep = prep
+  )
+  log_lik_weight(out, i = i, prep = prep)
+}
+
+log_lik_mixcure_loglogistic <- function(i, prep) {
+  shape <- get_dpar(prep, "shape", i = i)
+  mu <- get_dpar(prep, "mu", i)
+  inc <- get_dpar(prep, "inc", i)
+  args <- list(shape = shape, scale = mu, inc = inc)
+  out <- log_lik_censor(
+    dist = "mixcure_loglogistic", args = args, i = i, prep = prep
+  )
+  out <- log_lik_truncate(
+    out, cdf = pmixcure_loglogistic, args = args, i = i, prep = prep
+  )
+  log_lik_weight(out, i = i, prep = prep)
+}
+
 log_lik_zero_inflated_poisson <- function(i, prep) {
   zi <- get_dpar(prep, "zi", i)
   lambda <- get_dpar(prep, "mu", i)