# Negative binomial for catch data with GAM

I work with landing data (kg or tons), which don't follow a normal distribution and are more often than not over-dispersed.

I'm trying to obtain a standardized index of abundance using GAM to model the landings as a function of year, month, and region and using the effort (number of vessels) as an offset. I'm currently using a negative binomial with a log-link and it actually perform reasonably well, but I was wondering if there would be a more appropriate way to approach the problem.

I know the negative binomial is for count data, but I don't fully understand the implications of using it with continuous data. Also, shall I take the log of the effort when using it as an offset, and if yes, why so?

Here an example:

set.seed(123)

landings=round(rnegbin(100,theta=3)+runif(50, 0, 0.8),1)
nvessels=round(landings*0.2+runif(50,0,5))

mydata <- data.frame(landings=landings, year=rep(2020:2021,each=50), month=c(1:12,2:11,1:11,6:12,1:10,2:11,1:11,1:10,6:12,1:12), region=c(rep("A",12),rep("B",10),rep("C",11),rep("D",7),rep("E",10),rep("A",10),rep("B",11),rep("C",10),rep("D",7),rep("E",12)),effort=nvessels)

mydata <- mydata %>% mutate(year_f=as.factor(year),region_f=as.factor(region))

mod <- gam(landings ~ year_f + region_f + s(month), offset(effort), data=mydata, family=nb(link="log"), method="REML")



Apologies if this is a duplicate, but couldn't find an appropriate answer to this specific questions.

Thanks for any help.

Val

• You should certainly add your offset in the log scale because your response also uses a log link, the offset $t$ turns your model $log(\mu)=X\beta$ into $log(\mu / t)=X\beta$ or equivalently $log(\mu)=X\beta + log(t)$ (because $log(a/b)=log(a)-log(b)$). That's assuming your implementation doesn't take the log for you but I've never encountered that. Commented Nov 8, 2023 at 12:43

For things like this, the Tweedie (family of) distribution(s) is often used. The Tweedie distribution is a family of distributions and contains as special cases the Poisson and gamma distributions (for the parameter values allowed in mgcv; the family is richer than this but the distributions where the Tweedie power parameter $$p$$ is in the interval 1–2 are most typically used in applied work.)

The Tweedie family is defined for non-negative reals. One of the justifications of the Tweedie is that the response is the sum of events each yielding a real. Think of daily rainfall as the response, where the daily total is 0 or the sum of $$\ge 1$$ rainfall events each contributing >0 mm precipitation.

The landing data can be thought of as biomass. You could catch nothing (zero fish) and get $$y_i = 0$$ biomass, or for each landing you could catch $$F$$ fish, each having a non-negative weight, which sum to the total landing (biomass) observed as $$y_i$$. What is not observed are the weights of the individual fish. Hence the response is the aggregate of those $$F$$ weights and can take on any real value between 0 and infinity (theoretically).

Technically, the Tweedie with $$p = 1$$ is the Poisson and $$p = 2$$ the gamma distributions, so there is a transition point at somewhere as $$p$$ tends to 1 where the distribution becomes discrete.

To use this in gam(), we can modify your model to the following and fixing the offset specification as per @PBulls' comment

mod <- gam(landings ~ year_f + region_f + s(month) + offset(log(effort)),
data = mydata,
method = "REML")


If your data are overdispersed compared to the estimated conditional distribution (test this using say the DHARMa pkg, which will work for the tw() family and GAMs fitted with {mgcv} now), then you might need to go more exotic and model one or more of the other parameters of the Tweedie, namely the power $$p$$ and/or the scale $$\phi$$, using family twlss(), in which case your model might be (assuming you have sufficient data!)

mod <- gam(list(landings ~ year_f + region_f + s(month) + offset(log(effort)),
~ region_f,
~ year_f),
data = mydata,
family = twlss(),
method = "REML")


where I have (randomly) decided that the power parameter $$p$$ varies by region and the scale parameter $$\phi$$ varies by year. There's nothing stopping you including the same covariates in all three linear predictors, but that will require a lot of data!

• Could you comment in your answer why one should use the Tweedie distribution instead of negative binomial for this use-case? They are both combinations of gamma and poisson distributions, so I guess this has to do with the outcome being continuous vs. discrete. But I'm not familiar with "landing data" so I think that explanation could be really useful. Commented Nov 8, 2023 at 16:12
• The NB is defined for the set of non-negative integers, while the Tweedie family is defined for non-negative reals. One of the justifications of the Tweedie is that the response is the sum of events each yielding a real. Think of daily rainfall as the response, where the daily total is the sum or 0 or more rainfall events. Commented Nov 8, 2023 at 19:31
• Here the landing data is best thought of as biomass. You could catch nothing and get 0 biomass, or you can catch F fish, each of those F fish having a nonnegative weight, that sun to the total landing observed. What is not observed is the weights of the individual fish. Hence the response is the aggregate of those F weights and can take on any real value between 0 - infinity (theoretically). Technically, the Tweedie with $p =1$ is the Poisson & $p=2$ the gamma distributions, so there is a transition point at somewhere close to $p=1$ where the distribution becomes more discrete. It’s an odd one Commented Nov 8, 2023 at 19:37
• Thank you Gavin! That is very useful and much clearer. How many zeros do you think a Tweedie can handle before entering in the realm of zero-inflated models?
– Val
Commented Nov 12, 2023 at 12:00