From 489fbfb45f3d540ed7c19d2434abff640ae52cd3 Mon Sep 17 00:00:00 2001
From: Jillian Dunic <>
Date: Thu, 28 Mar 2024 11:16:32 -0700
Subject: [PATCH 1/6] add randomised quantile residuals for gengamma

---
 R/residuals.R                     | 24 ++++++++++++++++++++++++
 tests/testthat/test-3-families.R  |  1 -
 tests/testthat/test-5-residuals.R |  9 +++++++++
 3 files changed, 33 insertions(+), 1 deletion(-)

diff --git a/R/residuals.R b/R/residuals.R
index a3c9f65be..54e64c1db 100644
--- a/R/residuals.R
+++ b/R/residuals.R
@@ -147,6 +147,29 @@ qres_beta <- function(object, y, mu, ...) {
   stats::qnorm(u)
 }
 
+# Modified from flexsurv::pgengamma
+pgengamma <- function(q, mean, sigma, .Q, lower.tail = TRUE, log.p = FALSE) {
+  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
+  log_scale = mu - (sigma * sqrt(k) * log(k))
+  shape = 1 / (sigma * sqrt(k))
+  y <- log(q)
+  w <- (y - log_scale) * shape
+  prob <- pgamma(exp(w), shape = k)
+  if (!lower.tail) prob <- 1 - prob
+  if (log.p) prob <- log(prob)
+  prob
+}
+
+qres_gengamma <- function(object, y, mu, ...) {
+  theta <- get_pars(object)
+  .Q <- theta$gengamma_Q
+  sigma <- exp(theta$ln_phi)
+  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 +375,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 {
diff --git a/tests/testthat/test-3-families.R b/tests/testthat/test-3-families.R
index 511edd01b..80cbae8e0 100644
--- a/tests/testthat/test-3-families.R
+++ b/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")
 })
 
 
diff --git a/tests/testthat/test-5-residuals.R b/tests/testthat/test-5-residuals.R
index 1352e5300..df4e546c5 100644
--- a/tests/testthat/test-5-residuals.R
+++ b/tests/testthat/test-5-residuals.R
@@ -176,6 +176,15 @@ test_that("randomized quantile residuals work,", {
   )
   expect_error(residuals(fit), 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)
 })
 
 test_that("residuals() works", {

From 76818d4ee2ef64f164f48c2d17e1f2672564bb32 Mon Sep 17 00:00:00 2001
From: Jillian Dunic <>
Date: Tue, 2 Apr 2024 15:50:03 -0700
Subject: [PATCH 2/6] use non orig pgengamma to fix for Q <0

---
 R/residuals.R | 32 ++++++++++++++++++++------------
 1 file changed, 20 insertions(+), 12 deletions(-)

diff --git a/R/residuals.R b/R/residuals.R
index 54e64c1db..5cbf33258 100644
--- a/R/residuals.R
+++ b/R/residuals.R
@@ -147,25 +147,33 @@ qres_beta <- function(object, y, mu, ...) {
   stats::qnorm(u)
 }
 
-# Modified from flexsurv::pgengamma
+# 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) {
-  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
-  log_scale = mu - (sigma * sqrt(k) * log(k))
-  shape = 1 / (sigma * sqrt(k))
-  y <- log(q)
-  w <- (y - log_scale) * shape
-  prob <- pgamma(exp(w), shape = k)
-  if (!lower.tail) prob <- 1 - prob
-  if (log.p) prob <- log(prob)
-  prob
+  if (.Q != 0) {
+    y <- log(q)
+    beta <- .Q / sigma
+    k <- .Q^-2
+    if (.Q > 0) { a <- 1 } else { a <- -1 }
+    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 # What is expnu? where does nu fit?
+
+    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 = y, 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)
 }

From b5729bb71df0d6a1c6c13249e5c629e0bdd56271 Mon Sep 17 00:00:00 2001
From: Jillian Dunic <>
Date: Wed, 3 Apr 2024 12:56:30 -0700
Subject: [PATCH 3/6] 'nu' is in Stacy paramerisation

---
 R/residuals.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/residuals.R b/R/residuals.R
index 5cbf33258..06534c17c 100644
--- a/R/residuals.R
+++ b/R/residuals.R
@@ -156,7 +156,7 @@ pgengamma <- function(q, mean, sigma, .Q, lower.tail = TRUE, log.p = FALSE) {
     if (.Q > 0) { a <- 1 } else { a <- -1 }
     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 # What is expnu? where does nu fit?
+    expnu <- exp(.Q * w) * k
 
     if (.Q > 0) {
       stats::pgamma(q = expnu, shape = k, rate = 1, lower.tail = lower.tail, log.p = log.p)

From 48ed65a76adff7c549f0c6b8de892d1a76a6d591 Mon Sep 17 00:00:00 2001
From: Jillian Dunic <>
Date: Wed, 3 Apr 2024 13:30:14 -0700
Subject: [PATCH 4/6] include flexsurv to copyrights

---
 inst/COPYRIGHTS | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/inst/COPYRIGHTS b/inst/COPYRIGHTS
index ade0e0344..f11db4a8c 100644
--- a/inst/COPYRIGHTS
+++ b/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

From b100583dfc7b23ab96d27de2cd103e7259a9b456 Mon Sep 17 00:00:00 2001
From: Jillian Dunic <>
Date: Wed, 3 Apr 2024 16:10:19 -0700
Subject: [PATCH 5/6] fix typos

---
 R/residuals.R | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/R/residuals.R b/R/residuals.R
index 06534c17c..c1fa77bb0 100644
--- a/R/residuals.R
+++ b/R/residuals.R
@@ -153,7 +153,6 @@ pgengamma <- function(q, mean, sigma, .Q, lower.tail = TRUE, log.p = FALSE) {
     y <- log(q)
     beta <- .Q / sigma
     k <- .Q^-2
-    if (.Q > 0) { a <- 1 } else { a <- -1 }
     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
@@ -165,7 +164,7 @@ pgengamma <- function(q, mean, sigma, .Q, lower.tail = TRUE, log.p = FALSE) {
     }
   } else {
     # use lnorm with correction for Q = 0
-    stats::plnorm(q = y, meanlog = log(mean) - ((sigma^2) / 2), sdlog = sigma)
+    stats::plnorm(q = q, meanlog = log(mean) - ((sigma^2) / 2), sdlog = sigma)
   }
 }
 

From f0c8f3e00847f5fc2f69d524eb28ed062236586f Mon Sep 17 00:00:00 2001
From: Jillian Dunic <>
Date: Wed, 3 Apr 2024 16:10:39 -0700
Subject: [PATCH 6/6] add tests for Q < 0 and Q = 0

---
 tests/testthat/test-5-residuals.R | 29 +++++++++++++++++++++++++++++
 1 file changed, 29 insertions(+)

diff --git a/tests/testthat/test-5-residuals.R b/tests/testthat/test-5-residuals.R
index 799b1ea6a..1959204b5 100644
--- a/tests/testthat/test-5-residuals.R
+++ b/tests/testthat/test-5-residuals.R
@@ -188,6 +188,29 @@ test_that("randomized quantile residuals work,", {
     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", {
@@ -408,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)
+})
\ No newline at end of file