GLM or beta regression, with proportion data, and many 0 and 1

Question

I am trying to see if some anthropic variables (e.g., PopdensityAvg) explain animals' distribution. My dependent variable is the area occupied (non_null_cells) over the total area of interest (PAR), which can be expressed both as a percentage (range_perc, continuous variable, range: 0-100; e.g., 87.9) and as a proportion (range_prop, ranging 0-1, e.g. 0.87). Each species has also a status (Establishment, discrete variable; i.e., 0/1), and I expect to have different influences of each independent variable, based on species' status. For instance, for a species with status 0, I expect no influence of Roads, while for a species with status 1, I expect an influence of Roads. The values of the independent variables are averaged and measured throughout the species distribution range. Some of these independent variables can have a negative value (e.g., PopgrowthAvg). The data looks like this:

              Binomial Establishment      PAR non_null_cells range_prop range_perc PopdensityAvg PopgrowthAvg  Railways    Roads
1    Ammotragus_lervia             1 10840225       10762611  0.9928402   99.28402      68.68320    1.8325691 0.5568174 2.265508
2 Ammotragus_lervia_NE             0 10699087       10621473  0.9927457   99.27457      69.58296    1.8520451 0.5631284 2.288850
3  Apodemus_sylvaticus             1  2334776        2158612  0.9245478   92.45478     139.97713    2.6652732 0.9798067 2.927115
4     Atelerix_algirus             1   438917         383734  0.8742746   87.42746      76.02714    2.4077785 0.6011552 2.380663
5            Axis_axis             1     2072           1970  0.9507722   95.07722      19.40114    0.1917508 0.0000000 2.354487
6         Axis_axis_NE             0    18804          18804  1.0000000  100.00000      97.61798    2.7962866 0.3156403 2.715702

I got the recommendation of beta regression (from betareg package in R), but I have many 0's and 1's in range_prop. If I perform the suggested transformation (from an answer here Dealing with 0,1 values in a beta regression , "if y also assumes the extremes 0 and 1, a useful transformation in practice is (y * (n−1) + 0.5) / n where n is the sample size", it transforms my 1's in 0.95, and I don't know how good is that for my model, as I am very interested in outputs from observations with a range_prop value of 1).

I switched to GLM (package stats, I don't know if there are better options) thinking it would be better. I think I should use a Gamma distribution, based on the response variable (that is > 0 and is continuous). I understand that as my response variable's upper limit is 1 (or 100, if we consider range_perc), a beta regression would likely perform better compared to a Gamma distribution (that, if I am correct, assumes just values > 0, without an upper limit), but I don't know how to deal with so many 0's and 1's, and I do not understand if I should transform my data to comply with a particular modelling algorithm or if I can safely switch to more generalized techniques (such as GLMs). I am wondering if my reasonings are correct, and also how to deal with negative values in the independent variables (such as the ones in PopgrowthAvg, not present in the sample data I show here).

Edit: In my data, I also have the total area of interest, from which I extracted range_perc and range_prop for each species in Binomial. I added it to the sample data. I am using R 4.0.3

Edit n.2: After the useful answers, I tried a GLM with binomial distribution. As in Dunn & Smyth book on GLMs,

Binomial responses may be specified in the glm() formula in one of three ways:

1. The response can be supplied as the observed proportions yi, when the sample sizes mi are supplied as the weights in the call to glm().
2. The response can be given as a two-column array, the columns giving the numbers of successes and failures respectively in each group of size mi. The prior weights weights do not need to be supplied (r computes the weights m as the sum of the number of successes and failures for each row).

I tried both:

1. First approach

col_glm <- cbind(myData$non_null_cells, ((myData$PAR)-(myData$non_null_cells))) 
#use a 2 column matrix as the response variable with the first column being the counts of 'successes' and the second column being the counts of 'failures' 
#in this case, successes are the cells occupied by the species range (non_null_cells) and failures is the total area (PAR - non_null_cells)

glm2.2 <- glm(col_glm ~ PopdensityAvg + PopgrowthAvg + Railways + Roads,
              family = binomial,
              data = myData)

summary(glm2.2) 

Call:
glm(formula = col_glm ~ PopdensityAvg + PopgrowthAvg + Railways + 
    Roads, family = binomial, data = myData)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2612.51      2.49    129.26    610.45   1722.15  

Coefficients:
                 Estimate  Std. Error z value            Pr(>|z|)    
(Intercept)   -1.05252214  0.00230722  -456.2 <0.0000000000000002 ***
PopdensityAvg -0.05356983  0.00002347 -2282.4 <0.0000000000000002 ***
PopgrowthAvg   0.31591664  0.00027212  1160.9 <0.0000000000000002 ***
Railways       7.23302791  0.00304914  2372.2 <0.0000000000000002 ***
Roads          0.80775145  0.00109461   737.9 <0.0000000000000002 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 60080855  on 94  degrees of freedom
Residual deviance: 52299051  on 90  degrees of freedom
AIC: 52299902

Number of Fisher Scoring iterations: 6

2. Second approach

glm2.3 <- glm(range_prop ~ PopdensityAvg + PopgrowthAvg + Railways + Roads,
              weights = myData$PAR,
              family = binomial,
              data = myData) 
#specify the response as a proportion between 0 and 1, then specify another column as the 'weight' that gives the total number that the proportion is from 

summary(glm2.3)

Call:
glm(formula = range_prop ~ PopdensityAvg + PopgrowthAvg + Railways + 
    Roads, family = binomial, data = myData, weights = myData$PAR)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2612.51      2.49    129.26    610.45   1722.15  

Coefficients:
                 Estimate  Std. Error z value            Pr(>|z|)    
(Intercept)   -1.05252214  0.00230722  -456.2 <0.0000000000000002 ***
PopdensityAvg -0.05356983  0.00002347 -2282.4 <0.0000000000000002 ***
PopgrowthAvg   0.31591664  0.00027212  1160.9 <0.0000000000000002 ***
Railways       7.23302791  0.00304914  2372.2 <0.0000000000000002 ***
Roads          0.80775145  0.00109461   737.9 <0.0000000000000002 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 60080855  on 94  degrees of freedom
Residual deviance: 52299051  on 90  degrees of freedom
AIC: 52299902

Number of Fisher Scoring iterations: 6

The outputs are the same. I tried, with my naive knowledge, to interpret the results. The Deviance Residuals seem to have a high median value and 1Q and 3Q to be substantially apart, which I guess is not good.
Null deviance and Residual deviance are very high as well. Again: I guess it's not good.

Answer 1

The issues with deviance and potential perfect separation aren't troubling by themselves. It's not clear, however, that the analysis you describe will address your questions adequately.

You seem to have set up geographic grids, with each cell evaluated for presence or absence of a species. In a binomial model each cell would count as a trial. You have over 10 million trials (PAR) in each of your first 2 cases alone. The deviance necessarily increases with the total number of trials, so the large deviance values aren't an a priori problem.

Although your "fitted probabilities numerically 0 or 1 occurred" might be a sign of perfect separation, you don't need to worry about that, either, on its own. With your large number of trials, you don't have the usual problems of poor convergence, extreme regression coefficients and wide confidence intervals (CI). Your model converged in 6 iterations and shows narrow CI consistent with a very large number of trials.

The problem is that your trials aren't independent. The large numbers of trials within each species aren't independent of each other and can't necessarily be compared against those for the other species. Observations within a species will incorporate any species-specific effects that aren't captured in your model, effects that will differ among the species.

One species-specific characteristic is its own population density. With a random Poisson population distribution and an average population $\mu$ per cell, there is an expected fraction $\exp(-\mu)$ of empty cells. So the baseline probability of an empty cell, under a null hypothesis of a random distribution over the grid, depends on population density. That's not taken into account in your model. As it stands, your model might just be for a proxy of species' population densities.

A mixed model with random effects for species, for example via glmer() in the lme4 package, might help in this regard. Random intercepts would allow for differences among species in baseline probabilities of cell occupancy, associated with species' population densities (and other differences, like the ability to detect the species). I'm a bit hesitant, however, as a standard random effects model imposes a Gaussian distribution on those random effects, and that might not be appropriate for your data with over half of species having cell occupancies over 96.5% but a quarter below 59% and some at less than 1%.

Perhaps that skew in baseline probabilities is associated with corresponding skew in your predictors of interest and thus won't end up being a problem. Visual inspection of your data informed by your knowledge of the subject matter could help inform how to proceed.

If there are data available on the overall population density of each species, it would be best to include that in the model directly. I suspect that there are established ways to do this in the ecology literature, but that's outside my area of expertise.

Once you have a better form of model, there are a couple of other potential improvements to consider.

First, your current model assumes that the log-odds of occupancy increases linearly with each of your continuous predictors. That might not be a good assumption in general, and would seem particularly poor with a skewed predictor like your PopdensityAvg. Consider flexible modeling of the continuous predictors, with regression splines or another form of generalized additive model.

Second, your model as presented doesn't seem to address your intent:

For instance, for a species with status 0, I expect no influence of Roads, while for a species with status 1, I expect an influence of Roads.

To do that you need a model that includes status on its own and in an interaction with the other predictors whose associations with outcome you expect to differ depending on status.

Answer 2

A "quick" update: I looked into GAMMs and I think they might be a good option for my data. After exploring them a bit (and understanding they are more complicated), I thought I'd first try with GLMMs, to get a better understanding of the mixed part first. Again, my response variable is:

col_glm <- cbind(myData$non_null_cells, ((myData$PAR)-(myData$non_null_cells))) 

Binomial and Establishment are now factors (with, respectively, 64 and 2 levels). I first tried a model with the species (Binomial) as a random effect:

glmm_species <- glmer(col_glm ~ 
                        Establishment +
                PopdensityAvg + 
                PopgrowthAvg + 
                Railways + 
                Roads +
                (1 | Binomial),
              data = myData,
              family = binomial,
              nAGQ = 0,
              control = glmerControl(optimizer="bobyqa",tolPwrss=1e-3,optCtrl=list(maxfun=2e5)))

I added control = glmerControl(optimizer="bobyqa",tolPwrss=1e-3,optCtrl=list(maxfun=2e5)) and nAGQ = 0 after receiving warning messages regarding failed convergence and large eigenvalue ratios (model nearly unidentifiable). I also tried to rescale my predictors (e.g., scale(myData$PopdensityAvg)) or center the mean (e.g., myData$PopdensityAvg - mean(myData$PopdensityAvg)) but this didn't solve the problem. So I added the control = and nAGQ = parts. This is the summary of glmm_species:

Generalized linear mixed model fit by maximum likelihood (Adaptive Gauss-Hermite Quadrature, nAGQ = 0) ['glmerMod']
 Family: binomial  ( logit )
Formula: col_glm ~ Establishment + PopdensityAvg + PopgrowthAvg + Railways +      Roads + (1 | Binomial)
   Data: myData
Control: glmerControl(optimizer = "bobyqa", tolPwrss = 0.001, optCtrl = list(maxfun = 200000))

      AIC       BIC    logLik  deviance  df.resid 
1212169.9 1212187.8 -606077.9 1212155.9        88 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1006.34   -24.88     0.00    19.97   379.95 

Random effects:
 Groups   Name        Variance Std.Dev.
 Binomial (Intercept) 15.79    3.974   
Number of obs: 95, groups:  Binomial, 64

Fixed effects:
                  Estimate  Std. Error  z value             Pr(>|z|)    
(Intercept)     2.94643168  0.49753189    5.922        0.00000000318 ***
Establishment1 -0.18042578  0.00070243 -256.860 < 0.0000000000000002 ***
PopdensityAvg  -0.00634655  0.00004349 -145.917 < 0.0000000000000002 ***
PopgrowthAvg    0.20843341  0.00111517  186.907 < 0.0000000000000002 ***
Railways       -2.42287462  0.00582613 -415.863 < 0.0000000000000002 ***
Roads           0.93458051  0.00333284  280.415 < 0.0000000000000002 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) Estbl1 PpdnsA PpgrwA Ralwys
Estblshmnt1 -0.002                            
PopdnstyAvg  0.006 -0.047                     
PopgrwthAvg -0.009  0.040 -0.680              
Railways     0.001 -0.102 -0.440  0.056       
Roads       -0.011  0.099 -0.639  0.495 -0.187

I imagine Establishment status cannot be both a fixed and a random effect - indeed, if I try it, I receive a message of boundary (singular) fit, and the output of isSingular is TRUE, so I didn't try the same model but with (1 | Establishment).

I then tried to add, as suggested by EdM, the interaction between the predictors and the Establishment status:

glmm_species_est_interaction <- glmer(col_glm ~ 
                     Establishment +
                       PopdensityAvg + 
                     PopgrowthAvg + 
                     Railways + 
                     Roads + 
                     (Establishment * PopdensityAvg) + 
                     (Establishment * PopgrowthAvg) +
                     (Establishment * Railways) +
                     (Establishment * Roads) +
                       (1 | Binomial),
                   data = myData,
                   family = binomial,
                   nAGQ = 0,
                   control = glmerControl(optimizer="bobyqa",tolPwrss=1e-3,optCtrl=list(maxfun=2e5)))

Summary:

Generalized linear mixed model fit by maximum likelihood (Adaptive Gauss-Hermite Quadrature, nAGQ = 0) ['glmerMod']
 Family: binomial  ( logit )
Formula: col_glm ~ Establishment + PopdensityAvg + PopgrowthAvg + Railways +  
    Roads + (Establishment * PopdensityAvg) + (Establishment *  
    PopgrowthAvg) + (Establishment * Railways) + (Establishment *      Roads) + (1 | Binomial)
   Data: myData
Control: glmerControl(optimizer = "bobyqa", tolPwrss = 0.001, optCtrl = list(maxfun = 200000))

      AIC       BIC    logLik  deviance  df.resid 
 662632.4  662660.5 -331305.2  662610.4        84 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
    -778       -5        0        9 18455578 

Random effects:
 Groups   Name        Variance Std.Dev.
 Binomial (Intercept) 11.65    3.413   
Number of obs: 95, groups:  Binomial, 64

Fixed effects:
                               Estimate Std. Error  z value             Pr(>|z|)    
(Intercept)                   1.4670433  0.4312426    3.402             0.000669 ***
Establishment1                4.8699235  0.0157346  309.504 < 0.0000000000000002 ***
PopdensityAvg                -0.0070984  0.0001197  -59.318 < 0.0000000000000002 ***
PopgrowthAvg                  0.3039196  0.0034845   87.221 < 0.0000000000000002 ***
Railways                     -1.9641992  0.0118910 -165.183 < 0.0000000000000002 ***
Roads                         2.1504846  0.0088561  242.826 < 0.0000000000000002 ***
Establishment1:PopdensityAvg -0.0464354  0.0001039 -446.874 < 0.0000000000000002 ***
Establishment1:PopgrowthAvg   0.2266168  0.0026760   84.685 < 0.0000000000000002 ***
Establishment1:Railways       6.5647116  0.0198642  330.480 < 0.0000000000000002 ***
Establishment1:Roads         -2.3143838  0.0078887 -293.378 < 0.0000000000000002 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
                 (Intr) Estbl1 PpdnsA PpgrwA Ralwys Roads  Estblshmnt1:PpdA Estblshmnt1:PpgA Estblshmnt1:Rl
Estblshmnt1      -0.025                                                                                    
PopdnstyAvg       0.031 -0.447                                                                             
PopgrwthAvg      -0.032  0.511 -0.834                                                                      
Railways         -0.002 -0.127 -0.338  0.029                                                               
Roads            -0.036  0.650 -0.807  0.750 -0.119                                                        
Estblshmnt1:PpdA  0.009 -0.032  0.238 -0.269  0.021 -0.337                                                 
Estblshmnt1:PpgA  0.026 -0.556  0.658 -0.871  0.083 -0.599 -0.057                                          
Estblshmnt1:Rl   -0.019  0.484 -0.483  0.548 -0.326  0.671 -0.748           -0.357                         
Estblshmnt1:Rd    0.020 -0.936  0.347 -0.353  0.284 -0.611  0.057            0.377           -0.580  

This last model seems to perform slightly better (I am using AIC and F-test to check models' performance). However, I am unsure where to go from here, and how to interpret this output.