Merge pull request #484 from spsanderson/development

Fixes #470

Author: spsanderson

NAMESPACE

@@ -133,6 +133,7 @@ export(util_pareto_stats_tbl)
 export(util_poisson_aic)
 export(util_poisson_param_estimate)
 export(util_poisson_stats_tbl)
+export(util_rztnbinom_aic)
 export(util_t_stats_tbl)
 export(util_triangular_param_estimate)
 export(util_triangular_stats_tbl)
@@ -142,6 +143,7 @@ export(util_uniform_stats_tbl)
 export(util_weibull_aic)
 export(util_weibull_param_estimate)
 export(util_weibull_stats_tbl)
+export(util_ztn_binomial_param_estimate)
 importFrom(data.table,.SD)
 importFrom(data.table,as.data.table)
 importFrom(data.table,melt)

NEWS.md

@@ -5,6 +5,7 @@ None
 
 ## New Features
 1. #468 - Add function `util_negative_binomial_aic()` to calculate the AIC for the negative binomial distribution.
+2. #470 - Add function `util_ztn_binomial_param_estimate()` and `util_rztnbinom_aic()` to estimate the parameters and calculate the AIC for the zero-truncated negative binomial distribution.
 
 ## Minor Improvements and Fixes
 1. Fix #468 - Update `util_negative_binomial_param_estimate()` to add the use of

R/est-param-ztn-binmoial.R

@@ -0,0 +1,122 @@
+#' Estimate Zero Truncated Negative Binomial Parameters
+#'
+#' @family Parameter Estimation
+#' @family Binomial
+#' @family Zero Truncated Negative Distribution
+#'
+#' @author Steven P. Sanderson II, MPH
+#'
+#' @details This function will attempt to estimate the zero truncated negative
+#' binomial size and prob parameters given some vector of values.
+#'
+#' @description The function will return a list output by default, and  if the parameter
+#' `.auto_gen_empirical` is set to `TRUE` then the empirical data given to the
+#' parameter `.x` will be run through the `tidy_empirical()` function and combined
+#' with the estimated negative binomial data.
+#'
+#' One method of estimating the parameters is done via:
+#' -  MLE via \code{\link[stats]{optim}} function.
+#'
+#' @param .x The vector of data to be passed to the function.
+#' @param .auto_gen_empirical This is a boolean value of TRUE/FALSE with default
+#' set to TRUE. This will automatically create the `tidy_empirical()` output
+#' for the `.x` parameter and use the `tidy_combine_distributions()`. The user
+#' can then plot out the data using `$combined_data_tbl` from the function output.
+#'
+#' @examples
+#' library(dplyr)
+#' library(ggplot2)
+#' library(actuar)
+#'
+#' x <- as.integer(mtcars$mpg)
+#' output <- util_ztn_binomial_param_estimate(x)
+#'
+#' output$parameter_tbl
+#'
+#' output$combined_data_tbl |>
+#'   tidy_combined_autoplot()
+#'
+#' set.seed(123)
+#' t <- rztnbinom(100, 10, .1)
+#' util_ztn_binomial_param_estimate(t)$parameter_tbl
+#'
+#' @return
+#' A tibble/list
+#'
+#' @export
+#'
+
+util_ztn_binomial_param_estimate <- function(.x, .auto_gen_empirical = TRUE) {
+
+  # Check if actuar library is installed
+  if (!requireNamespace("actuar", quietly = TRUE)) {
+    stop("The 'actuar' package is needed for this function. Please install it with: install.packages('actuar')")
+  }
+
+  # Tidyeval ----
+  x_term <- as.numeric(.x)
+  sum_x <- sum(x_term, na.rm = TRUE)
+  minx <- min(x_term)
+  maxx <- max(x_term)
+  m <- mean(x_term, na.rm = TRUE)
+  n <- length(x_term)
+
+  # Negative log-likelihood function for optimization
+  nll_func <- function(params) {
+    size <- params[1]
+    prob <- params[2]
+    -sum(actuar::dztnbinom(x_term, size = size, prob = prob, log = TRUE))
+  }
+
+  # Initial parameter guesses
+  initial_params <- c(size = 1, prob = 0.5)  # Adjust based on your data
+
+  # Optimization using optim()
+  optim_result <- optim(initial_params, nll_func) |>
+    suppressWarnings()
+
+  # Extract estimated parameters
+  mle_size <- optim_result$par[1]
+  mle_prob <- optim_result$par[2]
+
+  # Create output tibble
+  ret <- tibble::tibble(
+    dist_type = "Zero-Truncated Negative Binomial",
+    samp_size = n,
+    min = minx,
+    max = maxx,
+    mean = m,
+    method = "MLE_Optim",
+    size = mle_size,
+    prob = mle_prob
+  )
+
+  # Attach attributes
+  attr(ret, "tibble_type") <- "parameter_estimation"
+  attr(ret, "family") <- "zero_truncated_negative_binomial"
+  attr(ret, "x_term") <- .x
+  attr(ret, "n") <- n
+
+  if (.auto_gen_empirical) {
+    # Generate empirical data
+    # Assuming tidy_empirical and tidy_combine_distributions functions exist
+    te <- tidy_empirical(.x = x_term)
+    td <- tidy_zero_truncated_negative_binomial(
+      .n = n,
+      .size = round(mle_size, 3),
+      .prob = round(mle_prob, 3)
+    )
+    combined_tbl <- tidy_combine_distributions(te, td)
+
+    output <- list(
+      combined_data_tbl = combined_tbl,
+      parameter_tbl = ret
+    )
+  } else {
+    output <- list(
+      parameter_tbl = ret
+    )
+  }
+
+  return(output)
+}

R/utils-aic-ztn-binomial.R

@@ -0,0 +1,93 @@
+#' Calculate Akaike Information Criterion (AIC) for Zero-Truncated Negative Binomial Distribution
+#'
+#' This function calculates the Akaike Information Criterion (AIC) for a
+#' zero-truncated negative binomial (ZTNB) distribution fitted to the provided data.
+#'
+#' @family Utility
+#' @author Steven P. Sanderson II, MPH
+#'
+#' @description
+#' This function estimates the parameters (`size` and `prob`) of a ZTNB
+#' distribution from the provided data using maximum likelihood estimation
+#' (via the `optim()` function), and then calculates the AIC value based on the
+#' fitted distribution.
+#'
+#' @param .x A numeric vector containing the data (non-zero counts) to be
+#'   fitted to a ZTNB distribution.
+#'
+#' @details
+#' **Initial parameter estimates:** The choice of initial values for `size`
+#' and `prob` can impact the convergence of the optimization. Consider using
+#' prior knowledge or method of moments estimates to obtain reasonable starting
+#' values.
+#'
+#' **Optimization method:** The default optimization method used is
+#' "Nelder-Mead". You might explore other optimization methods available in
+#' `optim()` for potentially better performance or different constraint
+#' requirements.
+#'
+#' **Data requirements:** The input data `.x` should consist of non-zero counts,
+#' as the ZTNB distribution does not include zero values.
+#'
+#' **Goodness-of-fit:** While AIC is a useful metric for model comparison, it's
+#' recommended to also assess the goodness-of-fit of the chosen ZTNB model using
+#' visualization (e.g., probability plots, histograms) and other statistical
+#' tests (e.g., chi-square goodness-of-fit test) to ensure it adequately
+#' describes the data.
+#'
+#' @examples
+#' library(actuar)
+#'
+#' # Example data
+#' set.seed(123)
+#' x <- rztnbinom(30, size = 2, prob = 0.4)
+#'
+#' # Calculate AIC
+#' util_rztnbinom_aic(x)
+#'
+#' @return The AIC value calculated based on the fitted ZTNB distribution to
+#'   the provided data.
+#'
+#' @name util_rztnbinom_aic
+NULL
+
+#' @export
+#' @rdname util_rztnbinom_aic
+
+util_rztnbinom_aic <- function(.x) {
+  # Check if actuar library is installed
+  if (!requireNamespace("actuar", quietly = TRUE)) {
+    stop("The 'actuar' package is needed for this function. Please install it with: install.packages('actuar')")
+  }
+
+  # Tidyeval
+  x <- as.numeric(.x)
+
+  # Get parameters
+  pe <- util_ztn_binomial_param_estimate(x)$parameter_tbl
+
+  # Negative log-likelihood function for zero-truncated negative binomial distribution
+  neg_log_lik_rztnbinom <- function(par, data) {
+    size <- par[1]
+    prob <- par[2]
+    -sum(actuar::dztnbinom(data, size = size, prob = prob, log = TRUE))
+  }
+
+  # Fit zero-truncated negative binomial distribution to data
+  fit_rztnbinom <- optim(
+    par = c(size = round(pe$size, 3), prob = round(pe$prob, 3)),
+    fn = neg_log_lik_rztnbinom,
+    data = x
+  ) |>
+    suppressWarnings()
+
+  # Extract log-likelihood and number of parameters
+  logLik_rztnbinom <- -fit_rztnbinom$value
+  k_rztnbinom <- 2  # Number of parameters (size and prob)
+
+  # Calculate AIC
+  AIC_rztnbinom <- 2 * k_rztnbinom - 2 * logLik_rztnbinom
+
+  # Return AIC value
+  return(AIC_rztnbinom)
+}