brms support for XBX regression

DESCRIPTION

@@ -2,8 +2,8 @@ Package: brms
 Encoding: UTF-8
 Type: Package
 Title: Bayesian Regression Models using 'Stan'
-Version: 2.22.3
-Date: 2024-10-30
+Version: 2.22.4
+Date: 2024-11-08
 Authors@R:
     c(person("Paul-Christian", "Bürkner", email = "paul.buerkner@gmail.com",
              role = c("aut", "cre")),
@@ -18,7 +18,8 @@ Authors@R:
       person("Hayden", "Rabel", role = c("ctb")),
       person("Simon C.", "Mills", role = c("ctb")),
       person("Stephen", "Wild", role = c("ctb")),
-      person("Ven", "Popov", role = c("ctb")))
+      person("Ven", "Popov", role = c("ctb")),
+      person("Ioannis", "Kosmidis", role = c("ctb")))
 Depends:
     R (>= 3.6.0),
     Rcpp (>= 0.12.0),
@@ -69,6 +70,7 @@ Suggests:
     statmod,
     digest,
     diffobj,
+    betareg,
     R.rsp,
     gtable,
     shiny,

NAMESPACE

@@ -627,6 +627,7 @@ export(von_mises)
 export(waic)
 export(weibull)
 export(wiener)
+export(xbeta)
 export(zero_inflated_beta)
 export(zero_inflated_beta_binomial)
 export(zero_inflated_binomial)

NEWS.md

@@ -1,3 +1,10 @@
+# brms 2.23
+
+### New Features
+
+* Fit extended-support Beta models via family `xbeta` 
+thanks to Ioannis Kosmidis. (#1698)
+
 # brms 2.22.0
 
 ### New Features

R/brm.R

@@ -4,7 +4,7 @@
 #' using Stan for full Bayesian inference. A wide range of distributions
 #' and link functions are supported, allowing users to fit -- among others --
 #' linear, robust linear, count data, survival, response times, ordinal,
-#' zero-inflated, hurdle, and even self-defined mixture models all in a
+#' zero-inflated, hurdle, extended-support beta regression, and even self-defined mixture models all in a
 #' multilevel context. Further modeling options include non-linear and
 #' smooth terms, auto-correlation structures, censored data, meta-analytic
 #' standard errors, and quite a few more. In addition, all parameters of the
@@ -412,7 +412,6 @@
 #' summary(fit7)
 #' conditional_effects(fit7)
 #'
-#'
 #' # use the future package for more flexible parallelization
 #' library(future)
 #' plan(multisession, workers = 4)

R/brmsformula.R

@@ -438,7 +438,7 @@
 #'   \code{Gamma}, \code{weibull}, \code{negbinomial}, and related zero-inflated
 #'   / hurdle families); \code{nu} (degrees of freedom parameter of the
 #'   \code{student} and \code{frechet} families); \code{phi} (precision
-#'   parameter of the \code{beta} and \code{zero_inflated_beta} families);
+#'   parameter of the \code{beta}, \code{zero_inflated_beta}, and \code{xbeta} families);
 #'   \code{kappa} (precision parameter of the \code{von_mises} family);
 #'   \code{beta} (mean parameter of the exponential component of the
 #'   \code{exgaussian} family); \code{quantile} (quantile parameter of the

R/distributions.R

@@ -1721,6 +1721,27 @@ pcox <- function(q, mu, bhaz, cbhaz, lower.tail = TRUE, log.p = FALSE) {
   out
 }
 
+# ensures that posterior_predict and friends can safely find these functions
+dxbeta <- function(...) {
+  require_package("betareg")
+  betareg::dxbeta(...)
+}
+
+pxbeta <- function(...) {
+  require_package("betareg")
+  betareg::pxbeta(...)
+}
+
+qxbeta <- function(...) {
+  require_package("betareg")
+  betareg::qxbeta(...)
+}
+
+rxbeta <- function(...) {
+  require_package("betareg")
+  betareg::rxbeta(...)
+}
+
 #' Zero-Inflated Distributions
 #'
 #' Density and distribution functions for zero-inflated distributions.

R/families.R

@@ -23,7 +23,7 @@
 #'   \code{hurdle_gamma}, \code{hurdle_lognormal}, \code{hurdle_cumulative},
 #'   \code{zero_inflated_binomial}, \code{zero_inflated_beta_binomial},
 #'   \code{zero_inflated_beta}, \code{zero_inflated_negbinomial},
-#'   \code{zero_inflated_poisson}, and \code{zero_one_inflated_beta}.
+#'   \code{zero_inflated_poisson}, \code{zero_one_inflated_beta}, and \code{xbeta}.
 #' @param link A specification for the model link function. This can be a
 #'   name/expression or character string. See the 'Details' section for more
 #'   information on link functions supported by each family.
@@ -94,6 +94,16 @@
 #'   \item{Families \code{beta}, \code{dirichlet}, and \code{logistic_normal}
 #'   can be used to model responses representing rates or probabilities.}
 #'
+#'   \item{Family \code{xbeta} extends the \code{beta} family to
+#'   support \code{[0, 1]} responses with exact \code{0}s and / or
+#'   \code{1}s, when each response takes values \code{0}, \code{1},
+#'   and \code{(0, 1)} according to a single process. If there is
+#'   merit in assuming that 0 and 1 values arise from different
+#'   processes than \code{(0, 1)} values, then the
+#'   \code{zero_inflated_beta}, \code{zero_one_inflated_beta} families
+#'   provide more flexibility.  For details see Kosmidis & Zeileis
+#'   (2024).}
+#'
 #'   \item{Family \code{asym_laplace} allows for quantile regression when fixing
 #'   the auxiliary \code{quantile} parameter to the quantile of interest.}
 #'
@@ -115,24 +125,26 @@
 #'   that cannot be explained by the primary distribution of the response.}
 #'   }
 #'
-#'   Below, we list all possible links for each family.
-#'   The first link mentioned for each family is the default.
-#'   \itemize{
-#'   \item{Families \code{gaussian}, \code{student}, \code{skew_normal},
-#'   \code{exgaussian}, \code{asym_laplace}, and \code{gen_extreme_value}
-#'   support the links (as names) \code{identity}, \code{log}, \code{inverse},
-#'   and \code{softplus}.}
+#'   Below, we list all possible links for each family.  The first
+#'   link mentioned for each family is the default.  \itemize{
+#'   \item{Families \code{gaussian}, \code{student},
+#'   \code{skew_normal}, \code{exgaussian}, \code{asym_laplace}, and
+#'   \code{gen_extreme_value} support the links (as names)
+#'   \code{identity}, \code{log}, \code{inverse}, and
+#'   \code{softplus}.}
 #'
 #'   \item{Families \code{poisson}, \code{negbinomial}, \code{geometric},
 #'   \code{zero_inflated_poisson}, \code{zero_inflated_negbinomial},
 #'   \code{hurdle_poisson}, and \code{hurdle_negbinomial} support
 #'   \code{log}, \code{identity}, \code{sqrt}, and \code{softplus}.}
 #'
-#'   \item{Families \code{binomial}, \code{bernoulli}, \code{beta_binomial},
-#'   \code{zero_inflated_binomial}, \code{zero_inflated_beta_binomial},
-#'   \code{Beta}, \code{zero_inflated_beta}, and \code{zero_one_inflated_beta}
-#'   support \code{logit}, \code{probit}, \code{probit_approx}, \code{cloglog},
-#'   \code{cauchit}, \code{identity}, and \code{log}.}
+#'   \item{Families \code{binomial}, \code{bernoulli},
+#'   \code{beta_binomial}, \code{zero_inflated_binomial},
+#'   \code{zero_inflated_beta_binomial}, \code{Beta},
+#'   \code{zero_inflated_beta}, \code{zero_one_inflated_beta}, and
+#'   \code{xbeta} support \code{logit}, \code{probit},
+#'   \code{probit_approx}, \code{cloglog}, \code{cauchit},
+#'   \code{identity}, and \code{log}.}
 #'
 #'   \item{Families \code{cumulative}, \code{cratio}, \code{sratio},
 #'   \code{acat}, and \code{hurdle_cumulative} support \code{logit},
@@ -174,7 +186,13 @@
 #'   have to use \code{brmsfamily} to specify the family with corresponding link
 #'   function.
 #'
-#' @seealso \code{\link[brms:brm]{brm}},
+#' @references
+#'
+#' Kosmidis I, Zeileis A (2024). Extended-Support Beta Regression for [0, 1] Responses.
+#' \emph{arXiv Preprint}. \doi{10.48550/arXiv.2409.07233}
+#'
+#' @seealso
+#'   \code{\link[brms:brm]{brm}},
 #'   \code{\link[stats:family]{family}},
 #'   \code{\link{customfamily}}
 #'
@@ -612,6 +630,15 @@ Beta <- function(link = "logit", link_phi = "log") {
               link_phi = link_phi)
 }
 
+#' @rdname brmsfamily
+#' @export
+xbeta <- function(link = "logit", link_phi = "log",
+                  link_kappa = "log") {
+  slink <- substitute(link)
+  .brmsfamily("xbeta", link = link, slink = slink,
+              link_phi = link_phi, link_kappa = link_kappa)
+}
+
 #' @rdname brmsfamily
 #' @export
 dirichlet <- function(link = "logit", link_phi = "log", refcat = NULL) {
@@ -1875,7 +1902,7 @@ family_bounds.brmsterms <- function(x, ...) {
     "hurdle_lognormal", "zero_inflated_poisson",
     "zero_inflated_negbinomial"
   )
-  beta_families <- c("beta", "zero_inflated_beta", "zero_one_inflated_beta")
+  beta_families <- c("beta", "zero_inflated_beta", "zero_one_inflated_beta", "xbeta")
   ordinal_families <- c("cumulative", "cratio", "sratio", "acat")
   if (family %in% pos_families) {
     out <- list(lb = 0, ub = Inf)

R/family-lists.R

@@ -61,6 +61,32 @@
   )
 }
 
+.family_xbeta <- function() {
+  list(
+    links = c(
+      "logit", "probit", "probit_approx", "cloglog",
+      "cauchit", "softit", "identity", "log"
+    ),
+    dpars = c("mu", "phi", "kappa"),
+    type = "real",
+    ybounds = c(0, 1),
+    closed = c(TRUE, TRUE),
+    ad = c("weights", "subset", "index"),
+    include = "fun_xbeta.stan",
+    prior = function(dpar, link = "identity", ...) {
+      if (dpar == "kappa") {
+        if (link == "identity") {
+          return("gamma(0.01, 0.01)")
+        } else {
+          return("student_t(3, 0, 2.5)")
+        }
+      }
+      NULL
+    },
+    normalized = ""
+  )
+}
+
 .family_beta_binomial <- function() {
   list(
     links = c(

R/log_lik.R

@@ -630,6 +630,19 @@ log_lik_beta <- function(i, prep) {
   log_lik_weight(out, i = i, prep = prep)
 }
 
+log_lik_xbeta <- function(i, prep) {
+  args <- list(
+    mu = get_dpar(prep, "mu", i = i),
+    phi = get_dpar(prep, "phi", i = i),
+    nu = get_dpar(prep, "kappa", i = i)
+  )
+  out <- log_lik_censor(dist = "xbeta", args = args, i = i, prep = prep)
+  out <- log_lik_truncate(
+    out, cdf = pxbeta, args = args, i = i, prep = prep
+  )
+  log_lik_weight(out, i = i, prep = prep)
+}
+
 log_lik_von_mises <- function(i, prep) {
   args <- list(
     mu = get_dpar(prep, "mu", i),

R/posterior_epred.R

@@ -450,6 +450,22 @@ posterior_epred_beta <- function(prep) {
   prep$dpars$mu
 }
 
+posterior_epred_xbeta <- function(prep) {
+  # see https://arxiv.org/abs/2409.07233 for details
+  mu <- get_dpar(prep, "mu")
+  phi <- get_dpar(prep, "phi")
+  nu <- get_dpar(prep, "kappa")
+  a <- mu * phi
+  b <- (1 - mu) * phi
+  d <- (1 + 2 * nu)
+  q0 <- nu / d
+  q1 <- (1 + nu) / d
+  t3 <- pbeta(q1, a, b)
+  t1 <- d * mu * (pbeta(q1, a + 1, b) - pbeta(q0, a + 1, b))
+  t2 <- nu * (t3 - pbeta(q0, a, b))
+  1 + t1 - t2 - t3
+}
+
 posterior_epred_von_mises <- function(prep) {
   prep$dpars$mu
 }

R/posterior_predict.R

@@ -682,6 +682,19 @@ posterior_predict_beta <- function(i, prep, ntrys = 5, ...) {
   )
 }
 
+posterior_predict_xbeta <- function(i, prep, ntrys = 5, ...) {
+  mu <- get_dpar(prep, "mu", i = i)
+  phi <- get_dpar(prep, "phi", i = i)
+  kappa <- get_dpar(prep, "kappa", i = i)
+  rcontinuous(
+    n = prep$ndraws,
+    dist = "xbeta",
+    mu = mu, phi = phi, nu = kappa,
+    lb = prep$data$lb[i], ub = prep$data$ub[i],
+    ntrys = ntrys
+  )
+}
+
 posterior_predict_von_mises <- function(i, prep, ntrys = 5, ...) {
   rcontinuous(
     n = prep$ndraws, dist = "von_mises",

R/priors.R

@@ -1119,7 +1119,7 @@ def_dpar_prior <- function(x, dpar) {
       sigma = def_scale_prior(x),
       shape = "student_t(3, 0, 2.5)",
       nu = "normal(2.7, 0.8)",
-      phi = "student_t(3, 0, 2.5)",
+      phi =  "student_t(3, 0, 2.5)",
       kappa = "normal(5.0, 0.8)",
       beta = "normal(1.7, 1.3)",
       zi = "logistic(0, 1)",

R/stan-likelihood.R

@@ -731,6 +731,12 @@ stan_log_lik_beta <- function(bterms, ...) {
   )
 }
 
+stan_log_lik_xbeta <- function(bterms, ...) {
+  reqn <- stan_log_lik_adj(bterms)
+  p <- stan_log_lik_dpars(bterms, reqn = reqn)
+  sdist("xbeta", p$mu, p$phi, p$kappa, vec = TRUE)
+}
+
 stan_log_lik_von_mises <- function(bterms, ...) {
   p <- stan_log_lik_dpars(bterms)
   sdist("von_mises", p$mu, p$kappa)