remove beta_binomial2 and replace with beta_binomial

open-AIMS · Sep 21, 2023 · c6c78fb · c6c78fb
1 parent bf8462d
commit c6c78fb
Show file tree

Hide file tree

Showing 38 changed files with 91 additions and 445 deletions.
diff --git a/NAMESPACE b/NAMESPACE
@@ -51,8 +51,6 @@ S3method(update,bnecfit)
 export(amend)
 export(average_estimates)
 export(bayesnecformula)
-export(beta_binomial2_lpmf)
-export(beta_binomial2_rng)
 export(bnec)
 export(bnec_newdata)
 export(bnf)
@@ -68,13 +66,10 @@ export(expand_manec)
 export(expand_nec)
 export(ggbnec_data)
 export(is_manecsummary)
-export(log_lik_beta_binomial2)
 export(make_brmsformula)
 export(models)
 export(nec)
 export(nsec)
-export(posterior_epred_beta_binomial2)
-export(posterior_predict_beta_binomial2)
 export(pull_brmsfit)
 export(pull_out)
 export(pull_prior)
@@ -86,6 +81,7 @@ importFrom(brms,as_draws_array)
 importFrom(brms,as_draws_df)
 importFrom(brms,bayes_R2)
 importFrom(brms,bernoulli)
+importFrom(brms,beta_binomial)
 importFrom(brms,brm)
 importFrom(brms,conditional_effects)
 importFrom(brms,fixef)
@@ -111,8 +107,6 @@ importFrom(dplyr,mutate)
 importFrom(dplyr,select)
 importFrom(evaluate,evaluate)
 importFrom(evaluate,is.warning)
-importFrom(extraDistr,dbbinom)
-importFrom(extraDistr,rbbinom)
 importFrom(formula.tools,`rhs<-`)
 importFrom(formula.tools,lhs)
 importFrom(formula.tools,rhs)

diff --git a/NEWS.md b/NEWS.md
@@ -91,7 +91,7 @@
 
 - Expanded `bnec`'s capacity to accept input x and y vectors, data frames ([`f256b39`](https://github.com/open-AIMS/bayesnec/commit/f256b399ab9115fffa7349a7a9daef21090f53f5)) and formulas ([`32e74ac`](https://github.com/open-AIMS/bayesnec/commit/32e74ac419c39c660aceb3d0914622de753a7a83)).
 
-- Corrected error for logit link cases for the beta_binomial2 where data contain 0 and 1 to ensure appropriate prior values on top and bot ([`4158237`](https://github.com/open-AIMS/bayesnec/commit/41582378a1a55c9420f69e578cfc98dc23182515)).
+- Corrected error for logit link cases for the "beta_binomial" where data contain 0 and 1 to ensure appropriate prior values on top and bot ([`4158237`](https://github.com/open-AIMS/bayesnec/commit/41582378a1a55c9420f69e578cfc98dc23182515)).
 
 - Series of internal fixes to standardise function class outputs ([`81369bb`](https://github.com/open-AIMS/bayesnec/commit/81369bbaef5e860410a5e2cc5227b6033687d36c), [`1c70efe`](https://github.com/open-AIMS/bayesnec/commit/1c70efeea54abe39c078ebfd014434e060c6f337), [`5bce2b5`](https://github.com/open-AIMS/bayesnec/commit/5bce2b5c40d8c1c480423529aaa59e0c82eda188), [`455ca70`](https://github.com/open-AIMS/bayesnec/commit/455ca70603a890b26a45b566975f21603f9f87df), [`5e6b41e`](https://github.com/open-AIMS/bayesnec/commit/5e6b41e6845321b5ff1f96c6733d59b6629fb707)).
 

diff --git a/R/amend.R b/R/amend.R
@@ -125,7 +125,7 @@ amend.bayesmanecfit <- function(object, drop, add, loo_controls, x_range = NA,
         x <- retrieve_var(bdat, "x_var", error = TRUE)
         y <- retrieve_var(bdat, "y_var", error = TRUE)
         custom_name <- check_custom_name(family)
-        if (family$family == "binomial" || custom_name == "beta_binomial2") {
+        if (family$family == "binomial" || family$family == "beta_binomial") {
           tr <- retrieve_var(bdat, "trials_var", error = TRUE)
           y <- y / tr
         }

diff --git a/R/autoplot.R b/R/autoplot.R
@@ -164,7 +164,7 @@ prep_raw_data <- function(brms_fit, bayesnecformula) {
   x_var <- attr(mod_dat, "bnec_pop")[["x_var"]]
   family <- brms_fit$family
   custom_name <- check_custom_name(family)
-  if (family$family == "binomial" | custom_name == "beta_binomial2") {
+  if (family$family == "binomial" | family == "beta_binomial") {
     trials_var <- attr(mod_dat, "bnec_pop")[["trials_var"]]
     r_df[[y_var]] <- r_df[[y_var]] / r_df[[trials_var]]
   } else {

diff --git a/R/bayesnec-package.R b/R/bayesnec-package.R
@@ -10,7 +10,7 @@
 #' @name bayesnec-package
 #' @aliases bayesnec
 #' @importFrom ggplot2 autoplot
-#' @importFrom brms bernoulli Beta negbinomial
+#' @importFrom brms bernoulli Beta negbinomial beta_binomial
 #'
 #' @references
 #' Bürkner P-C (2018) Advanced Bayesian Multilevel Modeling with the R Package

diff --git a/R/beta_binomial2.R b/R/beta_binomial2.R
diff --git a/R/bnec.R b/R/bnec.R
@@ -21,9 +21,7 @@
 #' ("fitting" and/or "weights"), each being a named \code{\link[base]{list}}
 #' containing the desired arguments to be passed on to \code{\link[brms]{loo}}
 #' (via "fitting") or to \code{\link[loo]{loo_model_weights}} (via "weights").
-#' If "fitting" is provided with argument \code{pointwise = TRUE}
-#' (due to memory issues) and \code{family = "beta_binomial2"}, the
-#' \code{\link{bnec}} will fail because that is a custom family. If "weights" is
+#' If "weights" is
 #' not provided by the user, \code{\link{bnec}} will set the default
 #' \code{method} argument in \code{\link[loo]{loo_model_weights}} to
 #' "pseudobma". See ?\code{\link[loo]{loo_model_weights}} for further info.
@@ -38,7 +36,7 @@
 #' to fit. See Details for more information.
 #' @param trials_var Removed in version 2.0. Use formula instead. Used to be a
 #' \code{\link[base]{character}} indicating the column
-#' heading for the number of "trials" for binomial or beta_binomial2 response
+#' heading for the number of "trials" for binomial or "beta_binomial" response
 #' data, as it appears in "data" (if data is supplied).
 #' @param random Removed in version 2.0. Use formula instead. Used to be a
 #' named \code{\link[base]{list}} containing the random model
@@ -124,8 +122,8 @@
 #' If not supplied via the \code{\link[brms]{brm}} argument \code{family}, the
 #' appropriate distribution will be guessed based on the characteristics of the
 #' input data. Guesses include: "bernoulli" / bernoulli / bernoulli(), "Beta" /
-#' Beta / Beta(), "binomial" / binomial / binomial(), "beta_binomial2" /
-#' beta_binomial2, "Gamma" / Gamma / Gamma(), "gaussian" / gaussian /
+#' Beta / Beta(), "binomial" / binomial / binomial(), "beta_binomial" /
+#' "beta_binomial", "Gamma" / Gamma / Gamma(), "gaussian" / gaussian /
 #' gaussian(), "negbinomial" / negbinomial / negbinomial(), or "poisson" /
 #' poisson / poisson(). Note, however, that "negbinomial" and "betabinomimal2"
 #' require knowledge on whether the data is over-dispersed. As

diff --git a/R/bnec_newdata.R b/R/bnec_newdata.R
@@ -53,7 +53,7 @@ bnec_newdata.bayesnecfit <- function(x, precision = 100, x_range = NA) {
   names(newdata) <- x_var
   fam_tag <- fit$family$family
   custom_name <- check_custom_name(fit$family)
-  if (fam_tag == "binomial" || custom_name == "beta_binomial2") {
+  if (fam_tag == "binomial" || fam_tag == "beta_binomial") {
     trials_var <- attr(data, "bnec_pop")[["trials_var"]]
     newdata[[trials_var]] <- 1
   }

diff --git a/R/check_data.R b/R/check_data.R
@@ -78,7 +78,7 @@ check_data <- function(data, family, model) {
     }
   }
   custom_name <- check_custom_name(family)
-  if (fam_tag == "binomial" || custom_name == "beta_binomial2") {
+  if (fam_tag == "binomial" || fam_tag == "beta_binomial") {
     mod_dat$trials <- retrieve_var(data, "trials_var", error = TRUE)
   }
   list(mod_dat = mod_dat, family = family)

diff --git a/R/check_models.R b/R/check_models.R
@@ -34,7 +34,7 @@ check_models <- function(model, family, data) {
     }
   }
   if (link_tag == "identity" & fam_tag %in%
-        c("bernoulli", "beta", "binomial", "custom")) {
+        c("bernoulli", "beta", "binomial", "beta_binomial")) {
     use_model <-  model[!model %in% c("neclin", "neclinhorme", "ecxlin")]
     drop_model <- setdiff(model, use_model)
     if (length(drop_model) > 0) {

diff --git a/R/data.R b/R/data.R
@@ -17,7 +17,7 @@ NULL
 #'
 #' @format An object of class \code{\link[brms]{customfamily}}
 #'
-#' @name beta_binomial2
+#' @name "beta_binomial"
 #' @docType data
 NULL
 

diff --git a/R/define_prior.R b/R/define_prior.R
@@ -18,17 +18,13 @@ define_prior <- function(model, family, predictor, response) {
   custom_name <- check_custom_name(family)
   if (link_tag %in% c("logit", "log")) {
     fam_tag <- "gaussian"
-  } else {
-    if (custom_name == "beta_binomial2") {
-      fam_tag <- custom_name
-    } else {
-      fam_tag <- family$family
-    }
-  }
-  if (custom_name == "beta_binomial2" || family$family == "binomial") {
+  } else { 
+    fam_tag <- family$family
+   }
+  if (family$family == "beta_binomial" || family$family == "binomial") {
     if (is.integer(response) || max(response) > 1) {
       stop("Response vector must be passed as a proportion to define_prior",
-           " (not as integers) for the binomial and beta_binomial2 families.")
+           " (not as integers) for the binomial and beta_binomial families.")
     }
   }
   response <- response_link_scale(response, family)
@@ -49,7 +45,7 @@ define_prior <- function(model, family, predictor, response) {
                                  ", ", sd(response) * 2.5, ")"),
                bernoulli = "beta(5, 2)",
                binomial = "beta(5, 2)",
-               beta_binomial2 = "beta(5, 2)",
+               "beta_binomial" = "beta(5, 2)",
                beta = "beta(5, 2)")
   y_b_prs <- c(Gamma = u_b_g,
                poisson = u_b_g,
@@ -59,7 +55,7 @@ define_prior <- function(model, family, predictor, response) {
                                  ", ", sd(response) * 2.5, ")"),
                bernoulli = "beta(2, 5)",
                binomial = "beta(2, 5)",
-               beta_binomial2 = "beta(2, 5)",
+               "beta_binomial" = "beta(2, 5)",
                beta = "beta(2, 5)")
   x_prs <- c(Beta = "beta(2, 2)",
              Gamma = paste0("gamma(5, ",
@@ -71,9 +67,9 @@ define_prior <- function(model, family, predictor, response) {
                                         probs = 0.5),
                                ", ", sd(predictor) * 10, ")"))
   lbs <- c(Gamma = 0, poisson = 0, negbinomial = 0, gaussian = NA,
-           bernoulli = 0, binomial = 0, beta_binomial2 = 0, beta = 0)
+           bernoulli = 0, binomial = 0, "beta_binomial" = 0, beta = 0)
   ubs <- c(Gamma = NA, poisson = NA, negbinomial = NA, gaussian = NA,
-           bernoulli = 1, binomial = 1, beta_binomial2 = 1, beta = 1)
+           bernoulli = 1, binomial = 1, "beta_binomial" = 1, beta = 1)
   # y-dependent priors
   pr_top <- prior_string(y_t_prs[fam_tag], nlpar = "top",
                          lb = lbs[fam_tag], ub = ubs[fam_tag])

diff --git a/R/expand_classes.R b/R/expand_classes.R
@@ -40,7 +40,7 @@ expand_nec <- function(object, formula, x_range = NA, precision = 1000,
   new_dat <- data.frame(x_seq)
   names(new_dat) <- x_var
   custom_name <- check_custom_name(fit$family)
-  if (fam_tag == "binomial" || custom_name == "beta_binomial2") {
+  if (fam_tag == "binomial" || fam_tag == "beta_binomial") {
     trials_col_name <- attr(mod_dat, "bnec_pop")[["trials_var"]]
     new_dat[[trials_col_name]] <- 1
   }

diff --git a/R/fit_bayesnec.R b/R/fit_bayesnec.R
@@ -42,14 +42,14 @@ fit_bayesnec <- function(formula, data, model = NA, brm_args,
       data[, y_var] <- y
       x_var <- bnec_pop_vars[[which(names(bnec_pop_vars) == "x_var")]]
       data[, x_var] <- x
-      if (family$family == "binomial" || custom_name == "beta_binomial2") {
+      if (family$family == "binomial" || family$family == "beta_binomial") {
         t_var <- bnec_pop_vars[[which(names(bnec_pop_vars) == "trials_var")]]
         data[, t_var] <- tr
       }
     }
   }
   custom_name <- check_custom_name(family)
-  if (family$family == "binomial" || custom_name == "beta_binomial2") {
+  if (family$family == "binomial" || family$family == "beta_binomial") {
     response <- y / tr
   } else {
     response <- y
@@ -58,15 +58,12 @@ fit_bayesnec <- function(formula, data, model = NA, brm_args,
   brm_args <- add_brm_defaults(brm_args, model, family, x, response,
                                skip_check, custom_name)
   all_args <- c(list(formula = brms_bf, data = quote(data)), brm_args)
-  if (custom_name == "beta_binomial2") {
-    all_args <- c(list(stanvars = stanvars), all_args)
-  }
   fit <- do.call(brm, all_args)
   pass <- are_chains_correct(fit, all_args$chains)
   if (!pass) {
     stop("Failed to fit model ", model, ".", call. = FALSE)
   }
-  msg_tag <- ifelse(family$family == "custom", custom_name, family$family)
+  msg_tag <- family$family
   message(paste0("Response variable modelled as a ", model, " model using a ",
                  msg_tag, " distribution."))
   out <- list(fit = fit, model = model, init = all_args$init,

diff --git a/R/helpers.R b/R/helpers.R
@@ -313,15 +313,7 @@ response_link_scale <- function(response, family) {
   lr <- linear_rescale
   custom_name <- check_custom_name(family)
   if (link_tag %in% c("logit", "log")) {
-    if (custom_name == "beta_binomial2") {
-      if (contains_zero(response)) {
-        response <- lr(response, r_out = c(min_z_val, max(response)))
-      }
-      if (contains_one(response)) {
-        response <- lr(response, r_out = c(min(response), max_o_val))
-      }
-      response <- binomial(link = link_tag)$linkfun(response)
-    } else if (family$family %in% c("bernoulli", "binomial")) {
+    if (family$family %in% c("bernoulli", "binomial", "beta_binomial")) {
       if (contains_zero(response)) {
         response <- lr(response, r_out = c(min_z_val, max(response)))
       }
@@ -484,7 +476,7 @@ add_brm_defaults <- function(brm_args, model, family, predictor, response,
     brm_args$prior <- priors
   }
   if (!("init" %in% names(brm_args)) || skip_check) {
-    msg_tag <- ifelse(family$family == "custom", custom_name, family$family)
+    msg_tag <- family$family
     message(paste0("Finding initial values which allow the response to be",
                    " fitted using a ", model, " model and a ", msg_tag,
                    " distribution."))

diff --git a/R/plot.R b/R/plot.R
@@ -94,7 +94,7 @@ plot.bayesnecfit <- function(x, ..., CI = TRUE, add_nec = TRUE,
   mod_dat <- model.frame(x$bayesnecformula, data = x$fit$data)
   y_var <- attr(mod_dat, "bnec_pop")[["y_var"]]
   x_var <- attr(mod_dat, "bnec_pop")[["x_var"]]
-  if (family == "binomial" | custom_name == "beta_binomial2") {
+  if (family == "binomial" | family == "beta_binomial") {
     trials_var <- attr(mod_dat, "bnec_pop")[["trials_var"]]
     y_dat <- x$fit$data[[y_var]] / x$fit$data[[trials_var]]
   } else {
@@ -250,7 +250,7 @@ plot.bayesmanecfit <- function(x, ..., CI = TRUE, add_nec = TRUE,
     x_var <- attr(bdat, "bnec_pop")[["x_var"]]
     family <- universal$fit$family$family
     custom_name <- check_custom_name(universal$fit$family)
-    if (family == "binomial" | custom_name == "beta_binomial2") {
+    if (family == "binomial" | family == "beta_binomial") {
       trials_var <- attr(bdat, "bnec_pop")[["trials_var"]]
       y_dat <- mod_dat[[y_var]] / mod_dat[[trials_var]]
     } else {

diff --git a/R/sysdata.rda b/R/sysdata.rda
diff --git a/R/validate_family.R b/R/validate_family.R
@@ -12,11 +12,7 @@ validate_family <- function(family) {
   if (inherits(family, "function")) {
     family <- family()
   } else if (is.character(family)) {
-    if (family == "beta_binomial2") {
-      family <- get(family)
-    } else {
       family <- get(family)(link = "identity")
-    }
   }
   if (!inherits(family, "family")) {
     stop("Argument \"family\" either is not an actual family, ",