Can GAM with Tweedie distribution be used for a binary response?

Question

I am trying species distribution models with GAM. The data set is the presence/absence of a species and environment variables in each grid. In the paper of Virgili et al. (2018), the authors suggested that, "to generate reliable habitat modelling predictions for rarely sighted species, using GAMs with a Tweedie distribution with presence-absence data". I tried the following code:

Model <- gam(Occurrence ~ s(Depth) + s(Temperature), data=dt, family=tw(link=log), method="REML")

The Model had a lower AIC and higher explained deviance than the model with binomial distribution. And the map of predict(Model, dt0, type="response") appeared fine.

However, as one can not opt for logit link function in Tweedie families, can I refer this predicted output as "probability of occurrence" as in the logistic regression?

Answer 1

It looks like they're actually modelling the abundance in that paper with NB, ZiP and TW distributions. They call it presence and absence but they really seem to be modelling densities as they say:

For each fitted model, predicted densities (in individuals per km2) were mapped on a 0.05°x0.05° resolution grid.

which would not be possible if they had reduced the data to 0s and 1s.

You shouldn't model binary (Bernoulli) data using a Tweedie. It might give a seemingly better predictive performance (using AIC) but there's nothing to stop it predicting a 0.5 or a 2 or any number greater than 1, which would be meaningless in terms of a binary response. So predictions from the model and uncertainty bands will likely be garbage unless you are very lucky indeed and happen to have covariate values that only predicted values in the range 0-1.

However, as one can not opt for logit link function in Tweedie families, can I refer this predicted output as "probability of occurrence" as in the logistic regression?

No, you cannot; the estimated response will be the expected rate of occurrence per unit effort (km² I think in the linked paper).

If you look at the output from qq.gam(Model, rep = 100) I'd be very surprised if your model could reproduce the data (i.e. the QQ plot of model residuals is contained within the reference band generated by sampling 100 new data sets for the response at the observed covariate values and residualising them.) I would expect the observed residuals to lie outside the bands generated by the model.

Another check you can do is via the DHARMa package. It will do something similar to what qq.gam() is doing but not quite the same; it generates randomised quantile residuals and then checks for deviations from expectations if the model were correct.

I should also add that it looks like these two data sets used in the paper were the result of a distance sampling survey and that means they should have considered the detectability of the individuals. This would require coupling a detection function to the GAM (which is very doable with mgcv: see Bravington et al; for a way to do this with mgcv), which they don't describe at all (so I presume they didn't do this).

Bravington, M.V., Miller, D.L., Hedley, S.L., 2021. Variance Propagation for Density Surface Models. J. Agric. Biol. Environ. Stat. https://doi.org/10.1007/s13253-021-00438-2