Zero- and one-inflated beta GAMM (Generalized additive mixed model) in mgcv

Question

I have vegetation cover (%) data [0,1] that includes 0's and 1's that I'd like to model with a beta GAMM, but don't understand the method for doing so. I've read that if the data includes 0's and 1's this makes it a "zero- and one- inflated model" which it involves adjusting eps within the family argument, but the instructions for what values that should take are not clear to me. All I can gather from other posts is that increasing eps within betar() needs to be done.

From the mgcv pdf: in betar(theta = NULL, link = "logit", eps=...), eps can be set to betar(theta = NULL, link = "logit",eps=.Machine$double.eps*100), but the context for this value is lost on me. Alternatively it can be set to eps <- 1e-10, which I assume is to ensure my values of 0 and 1 don't actually hit the [0,1] boundaries, but this doesn't make sense biologically. Maybe it doesn't need to?

I have many random effects that I need to include as this is a longitudinal study with fixed sites, within seasons, and years. My data has 903 rows, but only 3 values = 1 and zeros "can" be present (each value of "cover" is an average of 10 replicate quadrats of vegetation cover, so in theory, the mean of all 10 could be 0 at some point).

Example data:

# (Mean of 10 reps per site) Vegetation cover (%)
x <- as.data.frame(sample(seq(0.0, 1.00, 0.05), 1000, replace = TRUE))

# Sunlight
x$v1 <- sample(0:3000, 1000, replace = TRUE)

# Salinity
x$v2 <- rnorm(1000, mean = 25, sd = 5)

# Depth
x$v3 <- sample(30:200, 1000, replace = TRUE)

x$site <- seq(1, 50, 1)
x$site <- as.factor(x$site)

x$year <- rep(2011:2020, each=100)
x$year <- as.factor(x$year)
x$year.std <- x$year - min(x$year)    

x$season <- rep(c("DRY", "WET"), each=50)
x$season <- as.factor(x$season)

names(x)[1] <- "cover"


mod <- bam(cover ~ 
             s(v1) + 
             s(v2) + 
             s(v3) + 
             year.std * season + # Continuous trend probably varies by season
             
             s(site, bs = "re") +           # Non-independence of obs.
             s(site, year, bs = "re") +     # Sites within years
             s(site, year.std, bs = "re") + # Long-term trend by site
             s(season, year, bs = "re"),    # Seasons within years
           
           data = x, 
           method = 'fREML',
           discrete = TRUE,
           family=betar(link="logit"), # family=betar(link="logit", eps = ?),
           control = list(trace = TRUE))

Setup complete. Calling fit
Deviance = 1492.26942559717 Iterations - 1
Deviance = -891.84534192793 Iterations - 2
Deviance = -2184.96585537795 Iterations - 3 # R stalls here

Answer 1

What eps is doing is converting all 0s in the data to 0 + eps and the 1s to 1 - eps such that all data values are in the support of the distribution with finite density (at 0 and 1 the density of the beta distribution is infinite and as such isn't much use for modelling). Hence to use the common or garden variety beta distribution for data we have to adjust any 0s and 1s to be slightly not 0 or 1. If you have a few observations that are 0 or 1 then this fudge probably won't make much difference, but if they represent a reasonable proportion of your observations then things could go wrong badly.

If you more than just a few have true 0s and/or 1s, you'll need to use something other than mgcv and brms would be my choice.

(With vegatation cover, the way it is often collected by field ecologists can result in >100% because of the height structure of the vegatation. If that's the case here then I don't think any beta distribution would be suitable.)