pbs-assess / sdmTMB

:earth_americas: An R package for spatial and spatiotemporal GLMMs with TMB
https://pbs-assess.github.io/sdmTMB/
182 stars 26 forks source link

Gengamma #286

Closed James-Thorson-NOAA closed 4 months ago

James-Thorson-NOAA commented 6 months ago

Adding generalized-gamma distribution.

See test below, where the gengamma distribution results in gengamma_Q close to 0, i.e., collapsing to the lognormal distribution, and AIC confirms that it fits with slightly less than 2 AIC higher than the lognormal.

library(sdmTMB)
pcod_2011 = subset( pcod_2011, density > 0 )

# Build a mesh to implement the SPDE approach:
mesh <- make_mesh(pcod_2011, c("X", "Y"), cutoff = 20)

#
fit1 <- sdmTMB(
 density ~ 1 + depth_scaled,
 data = pcod_2011,
 mesh = mesh,
 spatial = "off",
 family = lognormal()
)
fit2 <- sdmTMB(
 density ~ 1 + depth_scaled,
 data = pcod_2011,
 mesh = mesh,
 spatial = "off",
 family = Gamma( link="log" )
)
fit3 <- sdmTMB(
 density ~ 1 + depth_scaled,
 data = pcod_2011,
 mesh = mesh,
 spatial = "off",
 family = gengamma()
)

AIC(fit1)
AIC(fit2)
AIC(fit3)
seananderson commented 6 months ago

Thanks!! Is it obvious what the distribution function would be? I.e. the equivalent of pgamma here so the quantile residuals work?

Also, I assume Q / lambda is the most intuitive parameter to report in print() rather than a derived version? As noted in the comments, the distribution approaches lognormal as Q -> 0.

Model fit by ML ['sdmTMB']
Formula: density ~ 1 + depth_scaled
Data: d
Family: gengamma(link = 'log')

             coef.est coef.se
(Intercept)      4.52    0.11
depth_scaled     0.20    0.11

Dispersion parameter: 1.41
Generalized gamma Q: 0.04
ML criterion at convergence: 2348.525
seananderson commented 6 months ago

Just noting for myself that there are some helpful examples here, which could be reparameterized for the Prentice version. Plus a bunch of implementations on CRAN to check against. First example https://search.r-project.org/CRAN/refmans/flexsurv/html/GenGamma.html

seananderson commented 6 months ago

Actually, I just realized it's the Prentice parameterization used in the flexsurv package already, so we can either import that pgengamma or use an open-source code snippet. https://search.r-project.org/CRAN/refmans/flexsurv/html/GenGamma.html