Merge pull request #333 from pbs-assess/add-gengamma-rqr

Add gengamma rqr

R/residuals.R

@@ -147,6 +147,36 @@ qres_beta <- function(object, y, mu, ...) {
   stats::qnorm(u)
 }
 
+# Modified from https://github.com/chjackson/flexsurv/blob/c765344fa798868036b841481fce2ea4d009d85e/src/gengamma.h#L63
+pgengamma <- function(q, mean, sigma, .Q, lower.tail = TRUE, log.p = FALSE) {
+  if (.Q != 0) {
+    y <- log(q)
+    beta <- .Q / sigma
+    k <- .Q^-2
+    mu <- log(mean) - lgamma((k * beta + 1) / beta) + lgamma(k) + log(k) / beta # Need to convert from mean to 'location' mu
+    w <- (y - mu) / sigma
+    expnu <- exp(.Q * w) * k
+
+    if (.Q > 0) {
+      stats::pgamma(q = expnu, shape = k, rate = 1, lower.tail = lower.tail, log.p = log.p)
+    } else {
+      stats::pgamma(q = expnu, shape = k, rate = 1, lower.tail = !lower.tail, log.p = log.p)
+    }
+  } else {
+    # use lnorm with correction for Q = 0
+    stats::plnorm(q = q, meanlog = log(mean) - ((sigma^2) / 2), sdlog = sigma)
+  }
+}
+
+qres_gengamma <- function(object, y, mu, ...) {
+  theta <- get_pars(object)
+  .Q <- theta$gengamma_Q
+  sigma <- exp(theta$ln_phi)
+  if (is_delta(object)) sigma <- sigma[2]
+  u <- pgengamma(q = y, mean = mu, sigma = sigma, .Q = .Q)
+  stats::qnorm(u)
+}
+
 #' Residuals method for sdmTMB models
 #'
 #' See the residual-checking vignette: `browseVignettes("sdmTMB")` or [on the
@@ -352,6 +382,7 @@ residuals.sdmTMB <- function(object,
       gamma_mix = qres_gamma_mix,
       lognormal_mix = qres_lognormal_mix,
       nbinom2_mix = qres_nbinom2_mix,
+      gengamma = qres_gengamma,
       cli_abort(paste(fam, "not yet supported."))
     )
   } else {

inst/COPYRIGHTS

@@ -5,6 +5,7 @@ The sdmTMB package includes other open source software components. The following
 is a list of these components.
 
 - VAST, https://github.com/James-Thorson-NOAA/VAST
+- flexsurv, https://github.com/chjackson/flexsurv
 - glmmTMB, https://github.com/glmmTMB/glmmTMB
 - SPDE barrier implementation, https://github.com/skaug/tmb-case-studies
 - mgcv, https://CRAN.R-project.org/package=mgcv
@@ -20,6 +21,9 @@ The following are licensed under GPL-3 as included for sdmTMB:
    bias correction (R/index.R get_generic()), and approach to optimization
    (modified into R/extra-optimization.R)
 
+2. flexsurv for the cumulative distribution function for the generalised
+   gamma distribution (modified into R/residuals.R pgengamma())
+
 The following are licensed under AGPL-3:
 
 1. glmmTMB for general approach to prediction (heavily modified into

tests/testthat/test-3-families.R

@@ -464,7 +464,6 @@ test_that("Generalized gamma works", {
   AIC(fit2)
   AIC(fit3)
 
-  expect_error(residuals(fit3), regexp = "supported")
 })
 
 

tests/testthat/test-5-residuals.R

@@ -176,8 +176,41 @@ test_that("randomized quantile residuals work,", {
     data = d, mesh = mesh,
     spatial = "off", spatiotemporal = "off"
   )
+
   expect_error(residuals(fit, type = "mle-mvn"), regexp = "truncated_nbinom1")
   check_resids_dharma(fit)
+
+  d <- sim_dat(gengamma())
+  fit <- sdmTMB(
+    observed ~ 1,
+    family = gengamma(),
+    data = d, mesh = mesh,
+    spatial = "off", spatiotemporal = "off"
+  )
+  check_resids(fit)
+
+  set.seed(1)
+  d <- sdmTMB_simulate(
+      formula = ~1,
+      data = predictor_dat,
+      mesh = mesh,
+      family = gengamma(),
+      range = 0.5,
+      phi = 0.2,
+      sigma_O = 0.2,
+      seed = 1,
+      B = 0,
+      control = sdmTMBcontrol(start = list(gengamma_Q = -1),
+                              map = list(gengamma_Q = factor(1))
+                              )
+  )
+  fit <- sdmTMB(
+      observed ~ 1,
+      family = gengamma(),
+      data = d, mesh = mesh,
+      spatial = "off", spatiotemporal = "off"
+  )
+  check_resids(fit)
 })
 
 test_that("residuals() works", {
@@ -398,3 +431,9 @@ test_that("old residual types get flagged with a message", {
   r <- residuals(fit, type = "mle-mvn")
   expect_length(r, 969L)
 })
+
+test_that("pgengamma works for Q = 0", {
+  set.seed(1)
+  p <- pgengamma(q = 1, mean = 0.5, sigma = 0.8, .Q = 0, lower.tail = TRUE, log.p = FALSE)
+  expect_equal(round(p, 4), 0.8973)
+})