Merge branch 'gengamma'

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: sdmTMB
 Title: Spatial and Spatiotemporal SPDE-Based GLMMs with 'TMB'
-Version: 0.4.3.9002
+Version: 0.4.2.9003
 Authors@R: c(
     person(c("Sean", "C."), "Anderson", , "[email protected]",
            role = c("aut", "cre"),

NAMESPACE

@@ -28,6 +28,7 @@ export(censored_poisson)
 export(delta_beta)
 export(delta_gamma)
 export(delta_gamma_mix)
+export(delta_gengamma)
 export(delta_lognormal)
 export(delta_lognormal_mix)
 export(delta_poisson_link_gamma)
@@ -38,6 +39,7 @@ export(dharma_residuals)
 export(extract_mcmc)
 export(gamma_mix)
 export(gather_sims)
+export(gengamma)
 export(get_cog)
 export(get_crs)
 export(get_index)

NEWS.md

@@ -1,5 +1,9 @@
 # sdmTMB (development version)
 
+* Add the generalized gamma distribution (thanks to J.T. Thorson with additional
+  work by J.C. Dunic.) See `gengamma()`. This distribution is still in a testing
+  phase and is not recommended for applied use yet. #286
+
 * Detect possible issue with factor(time) in formula if same column name is used
   for `time` and `extra_time` is specified. #320
 

R/enum.R

@@ -16,7 +16,8 @@
   censored_poisson  = 12,
   gamma_mix = 13,
   lognormal_mix = 14,
-  nbinom2_mix = 15
+  nbinom2_mix = 15,
+  gengamma = 16
 )
 .valid_link <- c(
   identity = 0,

R/families.R

@@ -72,6 +72,40 @@ lognormal <- function(link = "log") {
   add_to_family(x)
 }
 
+#' @export
+#' @rdname families
+#' @examples
+#' gengamma(link = "log")
+#' @details
+#' The `gengamma()` family was implemented by J.T. Thorson and uses the Prentice
+#' (1974) parameterization such that the lognormal occurs as the internal
+#' parameter `gengamma_Q` (reported in `print()` as "Generalized gamma lambda")
+#' approaches 0.
+#'
+#' @references
+#'
+#' *Generalized gamma family*:
+#'
+#' Prentice, R.L. 1974. A log gamma model and its maximum likelihood estimation.
+#' Biometrika 61(3): 539–544. \doi{10.1093/biomet/61.3.539}
+#'
+#' Stacy, E.W. 1962. A Generalization of the Gamma Distribution. The Annals of
+#' Mathematical Statistics 33(3): 1187–1192. Institute of Mathematical
+#' Statistics.
+
+gengamma <- function(link = "log") {
+  linktemp <- substitute(link)
+  if (!is.character(linktemp))
+    linktemp <- deparse(linktemp)
+  okLinks <- c("identity", "log", "inverse")
+  if (linktemp %in% okLinks)
+    stats <- stats::make.link(linktemp)
+  else if (is.character(link))
+    stats <- stats::make.link(link)
+  x <- c(list(family = "gengamma", link = linktemp), stats)
+  add_to_family(x)
+}
+
 #' @details The families ending in `_mix()` are 2-component mixtures where each
 #'   distribution has its own mean but a shared scale parameter.
 #'   (Thorson et al. 2011). See the model-description vignette for details.
@@ -328,6 +362,30 @@ delta_gamma_mix <- function(link1 = "logit", link2 = "log") {
        clean_name = "delta_gamma_mix(link1 = 'logit', link2 = 'log')"), class = "family")
 }
 
+#' @export
+#' @examples
+#' delta_gengamma()
+#' @rdname families
+delta_gengamma <- function(link1 = "logit", link2 = "log", type = c("standard", "poisson-link")) {
+  type <- match.arg(type)
+  l1 <- substitute(link1)
+  if (!is.character(l1)) l1 <- deparse(l1)
+  l2 <- substitute(link2)
+  if (!is.character(l2)) l2 <- deparse(l2)
+  f1 <- binomial(link = l1)
+  f2 <- gengamma(link = l2)
+  if (type == "poisson-link") {
+    .type <- "poisson_link_delta"
+    clean_name <- paste0("delta_gengamma(link1 = '", l1, "', link2 = '", l2, "', type = 'poisson-link')")
+  } else {
+    .type <- "standard"
+    clean_name <- paste0("delta_gengamma(link1 = '", l1, "', link2 = '", l2, "')")
+  }
+  structure(list(f1, f2, delta = TRUE, link = c("logit", "log"),
+    type = .type, family = c("binomial", "gengamma"),
+    clean_name = clean_name), class = "family")
+}
+
 #' @export
 #' @examples
 #' delta_lognormal()

R/fit.R

@@ -27,8 +27,10 @@ NULL
 #'   \code{\link[sdmTMB:families]{censored_poisson()}},
 #'   \code{\link[sdmTMB:families]{gamma_mix()}},
 #'   \code{\link[sdmTMB:families]{lognormal_mix()}},
-#'   \code{\link[sdmTMB:families]{student()}}, and
-#'   \code{\link[sdmTMB:families]{tweedie()}}. Supports the delta/hurdle models:
+#'   \code{\link[sdmTMB:families]{student()}},
+#'   \code{\link[sdmTMB:families]{tweedie()}}, and
+#'   \code{\link[sdmTMB:families]{gengamma()}}.
+#'   Supports the delta/hurdle models:
 #'   \code{\link[sdmTMB:families]{delta_beta()}},
 #'   \code{\link[sdmTMB:families]{delta_gamma()}},
 #'   \code{\link[sdmTMB:families]{delta_gamma_mix()}},
@@ -1151,6 +1153,7 @@ sdmTMB <- function(
     ln_kappa   = matrix(0, 2L, n_m),
     # ln_kappa   = rep(log(sqrt(8) / median(stats::dist(spde$mesh$loc))), 2),
     thetaf     = 0,
+    gengamma_Q = 1, # Not defined at exactly 0
     logit_p_mix = 0,
     log_ratio_mix = 0,
     ln_phi     = rep(0, n_m),
@@ -1182,6 +1185,7 @@ sdmTMB <- function(
   tmb_map$b_j <- NULL
   if (delta) tmb_map$b_j2 <- NULL
   if (family$family[[1]] == "tweedie") tmb_map$thetaf <- NULL
+  if ("gengamma" %in% family$family) tmb_map$gengamma_Q <- NULL
   if (family$family[[1]] %in% c("gamma_mix", "lognormal_mix", "nbinom2_mix")) {
     tmb_map$log_ratio_mix <- NULL
     tmb_map$logit_p_mix <- NULL

R/print.R

@@ -261,9 +261,13 @@ print_other_parameters <- function(x, m = 1L) {
   b <- tidy(x, "ran_pars", model = m, silent = TRUE)
 
   get_term_text <- function(term_name = "", pretext = "") {
+    b2 <- as.list(x$sd_report, what = "Estimate")
     if (term_name %in% b$term) {
       a <- mround(b$estimate[b$term == term_name], 2L)
       a <- paste0(pretext, ": ", a, "\n")
+    } else if (term_name %in% names(b2)) {
+      a <- mround(b2[[term_name]], 2L)
+      a <- paste0(pretext, ": ", a, "\n")
     } else {
       a <- ""
     }
@@ -272,6 +276,9 @@ print_other_parameters <- function(x, m = 1L) {
 
   phi <- get_term_text("phi", "Dispersion parameter")
   tweedie_p <- get_term_text("tweedie_p", "Tweedie p")
+  gengamma_par <- if ('gengamma' %in% family(x)[[m]]) {
+    get_term_text("gengamma_Q", "Generalized gamma lambda") 
+    } else ""
   sigma_O <- get_term_text("sigma_O", "Spatial SD")
   xtra <- if (x$spatiotemporal[m] == "ar1") "marginal " else ""
   sigma_E <- get_term_text("sigma_E",
@@ -292,7 +299,7 @@ print_other_parameters <- function(x, m = 1L) {
     sigma_Z <- ""
   }
 
-  named_list(phi, tweedie_p, sigma_O, sigma_E, sigma_Z, rho)
+  named_list(phi, tweedie_p, sigma_O, sigma_E, sigma_Z, rho, gengamma_par)
 }
 
 print_header <- function(x) {
@@ -340,6 +347,7 @@ print_one_model <- function(x, m = 1) {
 
   cat(other$phi)
   cat(other$tweedie_p)
+  cat(other$gengamma_par)
   cat(other$rho)
   cat(range)
   cat(other$sigma_O)

src/sdmTMB.cpp

@@ -20,7 +20,8 @@ enum valid_family {
   censored_poisson_family  = 12,
   gamma_mix_family = 13,
   lognormal_mix_family = 14,
-  nbinom2_mix_family = 15
+  nbinom2_mix_family = 15,
+  gengamma_family = 16
 };
 
 enum valid_link {
@@ -233,6 +234,7 @@ Type objective_function<Type>::operator()()
   PARAMETER_ARRAY(ln_kappa);    // Matern parameter
 
   PARAMETER(thetaf);           // tweedie only
+  PARAMETER(gengamma_Q);           // gengamma only
   PARAMETER(logit_p_mix);           // ECE / positive mixture only
   PARAMETER(log_ratio_mix);           // ECE / positive mixture only
 
@@ -934,6 +936,11 @@ Type objective_function<Type>::operator()()
           }
           break;
         }
+          case gengamma_family: {
+            tmp_ll = sdmTMB::dgengamma(y_i(i,m), mu_i(i,m), phi(m), gengamma_Q, true);
+            SIMULATE{y_i(i,m) = sdmTMB::rgengamma(mu_i(i,m), phi(m), gengamma_Q);}
+            break;
+          }
         default:
           error("Family not implemented.");
         }

src/utils.h

@@ -4,6 +4,61 @@ bool isNA(Type x) {
   return R_IsNA(asDouble(x));
 }
 
+// dgengamma
+// Written by J. Thorson based on scripts listed below
+// using Prentice-1974 parameterization for lambda instead of k, so that lognormal occurs as lambda -> 0
+// using mean parameterization to back out theta
+// CV is a function of sigma and lambda and NOT mean (i.e., CV is fixed for all values of mean)
+// See: C:\Users\James.Thorson\Desktop\Work files\AFSC\2021-10 -- Generalized gamma-lognormal\Explore gengamma.R
+template<class Type>
+Type dgengamma( Type x,
+                Type mean,
+                Type sigma,
+                Type Q,
+                int give_log=0){
+
+
+  // First convert from mean to mu
+  // Rename based on flexdist
+  Type lambda = Q;
+  // Using https://stats.stackexchange.com/questions/345564/generalized-gamma-log-normal-as-limiting-special-case
+  Type k = pow( lambda, -2);
+  Type beta = pow( sigma, -1) * lambda;
+  // Use wikipedia expression for the mean:   mean = a * gamma((d+1)/p) / gamma(d/p) where d/p = k, theta = a, and beta = p
+  // https://en.wikipedia.org/wiki/Generalized_gamma_distribution#Software_implementation
+  Type log_theta = log(mean) - lgamma( (k*beta+1)/beta ) + lgamma( k );
+  // Using https://stats.stackexchange.com/questions/345564/generalized-gamma-log-normal-as-limiting-special-case
+  Type mu = log_theta + log(k) / beta;
+
+  // Next evaluate PDF
+  // from https://github.com/chjackson/flexsurv-dev/blob/master/src/gengamma.h#L54-L56
+  Type y = log(x);
+  Type w = (y - mu) / sigma;
+  Type qi = pow(Q, -2);
+  Type qw = Q * w;                 // 0.5*log(pow(x,2)) as trick for abs(log(x))
+  Type logres = -log(sigma*x) + 0.5*log(pow(lambda,2)) * (1 - 2 * qi) + qi * (qw - exp(qw)) - lgamma(qi);
+
+  // return stuff
+  if(give_log) return logres; else return exp(logres);
+}
+
+// rgengamma
+// Written by J. Thorson based on scripts listed below
+template<class Type>
+Type rgengamma( Type mean,
+                Type sigma,
+                Type Q){
+
+  // See: C:\Users\James.Thorson\Desktop\Work files\AFSC\2021-10 -- Generalized gamma-lognormal\Explore gengamma.R
+  Type lambda = Q;
+  Type k = pow( lambda, -2 );
+  Type beta = pow( sigma, -1 ) * lambda;
+  Type log_theta = log(mean) - lgamma( (k*beta+1)/beta ) + lgamma( k );
+  Type w = log(rgamma(k, Type(1.0)));
+  Type y = w/beta + log_theta;
+  return exp(y);
+}
+
 template <class Type>
 Type ppois_log(Type x, Type lambda) {
   return atomic::Rmath::Rf_ppois(asDouble(x), asDouble(lambda), true, true);