Small Sample Sizes and Zero Inflated Count Data in R

Question

I am working to produce a model in R for seed germination count data with lots of zeros (around 50% of the 264 total observations). The purpose is to determine the effect treatments have on plant species as a whole and individually. Originally, I had fit a generalised linear mixed effects model (binomial error structure) and though I got decent results for an entire dataset model, when I looked at individual species the model was unable to capture a seemingly obvious effect of a treatment when all the alternative cases were 0. For example the below is the summary output of the model when only filtered to one seed with around 24 observations. The sample data behind the output is below as well.

n<- c(19,20,20,20,20,20,19,19,20,20,20,20,19,20,20,20,20,20,20,20,19,19,19,20)
count<-c(8,11,13,13,15,11,0,0,0,0,0,0,9,9,13,11,14,8,0,0,0,0,0,0)
GA<-c(1,1,1,1,1,1,0,0,0,0,0,0,1,1,1,1,1,1,0,0,0,0,0,0)
Sm<-c(1,1,1,0,0,0,1,1,1,0,0,0,1,1,1,0,0,0,1,1,1,0,0,0)
light<-c(1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0)
mod_data<-data.frame(n, count, GA, Sm, light)

Call:
glm(formula = cbind(count, n - count) ~ GA + Sm + light, 
    family = binomial, data = mod_data)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)  -23.7210  5592.9950  -0.004    0.997
GA            24.0074  5592.9950   0.004    0.997
Sm            -0.2706     0.2627  -1.030    0.303
light          0.2410     0.2627   0.917    0.359

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 251.9059  on 23  degrees of freedom
Residual deviance:   9.9247  on 20  degrees of freedom
AIC: 58.669

Number of Fisher Scoring iterations: 20 

I originally expected this seed to have a significant GA treatment as looking at the raw results, all the GA observations produced higher counts than their non-GA results (of all zero). However, I have since discovered that this type of model doesn't utilise zeros and therefore cannot produce a significant effect of GA (even though the overall model does).

This led me to investigate zero-inflated models like the zero-inflated negative binomial GLM (using the pscl package). However, when I looked at a similar model for that same seed I received an warning message and all NAs for my treatments:

Warning message:
In value[[3L]](cond) :
  system is computationally singular: reciprocal condition number = 1.09697e-23FALSE
> summary(mod)

Call:
zeroinfl(formula = count ~ GA + Sm + light | GA + Sm + light, data = mod_data, dist = "negbin")

Pearson residuals:
       Min         1Q     Median         3Q        Max 
-1.001e+00 -3.804e-02 -6.695e-11 -6.263e-11  9.649e-01 

Count model coefficients (negbin with log link):
            Estimate Std. Error z value Pr(>|z|)
(Intercept) -21.2880         NA      NA       NA
GA           23.7197         NA      NA       NA
Sm           -0.1335         NA      NA       NA
light         0.1038         NA      NA       NA
Log(theta)   15.5141         NA      NA       NA

Zero-inflation model coefficients (binomial with logit link):
              Estimate Std. Error z value Pr(>|z|)
(Intercept)  2.557e+01         NA      NA       NA
GA          -5.113e+01         NA      NA       NA
Sm          -9.764e-09         NA      NA       NA
light       -9.764e-09         NA      NA       NA

Theta = 5466023.2431 
Number of iterations in BFGS optimization: 20 
Log-likelihood: -27.84 on 9 Df

Its not entirely clear to me why NAs are being produced here however, the overall model is again producing decent (and better than first model) results.

My question is, is there a better approach in examining the effect of treatments on an individual seed like this? Or is it better to just look at the overall models and initial data exploration for the individual seeds? Thanks in advance.

Answer 1

The problem you are observing with lack of significance of GA has nothing to do with zero inflation or with random effects. It is simply a limitation of Wald tests for count models. If you replace the Wald tests with likelihood ratio tests, then the problem disappears:

> n <-  c(19,20,20,20,20,20,19,19,20,20,20,20,19,20,20,20,20,20,20,20,19,19,19,20)
> count <- c(8,11,13,13,15,11,0,0,0,0,0,0,9,9,13,11,14,8,0,0,0,0,0,0)
> GA <- c(1,1,1,1,1,1,0,0,0,0,0,0,1,1,1,1,1,1,0,0,0,0,0,0)
> Sm <- c(1,1,1,0,0,0,1,1,1,0,0,0,1,1,1,0,0,0,1,1,1,0,0,0)
> light <- c(1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0)
> 
> p <- count/n
> fit <- glm(p~GA+Sm+light,family=binomial(),weights=n)
> anova(fit,test="Chisq")
Analysis of Deviance Table

Model: binomial, link: logit

Response: p

Terms added sequentially (first to last)


      Df Deviance Resid. Df Resid. Dev Pr(>Chi)    
NULL                     23    251.906             
GA     1  240.079        22     11.827   <2e-16 ***
Sm     1    1.060        21     10.768   0.3033    
light  1    0.843        20      9.925   0.3585    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The problem here is that the GA variable perfectly separates the zero count responses from the non-zero. Hence the maximum likelihood estimate of the coefficient for GA in the linear model is infinite and so is its standard error. The z-statistic, which is the ratio of the coefficient to its standard error, is therefore undefined. In floating point arithmetic, the coefficient never becomes infinite, but the standard error approaches infinity faster than the coefficient as the iteration proceeds, so the z-statistic converges to zero and appears to be non-significant.

The solution is simply not to use z-statistics. The likelihood ratio tests depend only on fitted values, which do not become infinite, and hence remain defined. For your data, the likelihood ratio test statistic for testing the GA coefficient equal to zero (minus twice the log-likelihood difference) is 240.079. Under the null hypothesis, this statistic follows a chisquare distribution on 1 df, so the observed value of 240.079 is astronomically significant.

Note that the other two covariates Sm and light have finite coefficients. For those covariates, the Wald test p-values and likelihood ratio p-values are almost the same.

The same problem occurs with Wald tests for binomial, Poisson and negative binomial generalized linear models with logit or log-link functions when one of the groups being compared has all zero counts and the corresponding coefficient becomes infinite. The problem has been observed many times over the years but strangely the solution never seems to have become well known. We discuss the problem and solution for binomial glms in Section 9.9 ("When Wald tests fail) of Dunn & Smyth (2018). We discuss the same problem for negative binomial models in Section 7.1 of Robinson & Smyth (2008). In that paper we show that the likelihood ratio test p-values are a reasonable approximation to exact tests but are slightly anticonservative. Rainey (2024) is another paper giving the same advice that I have always given. Greenland (1992) also recommends likelihood-ratio tests for logistic regression when the counts are small.

Reference

Dunn PK, Smyth GK (2018). Generalized linear models with examples in R. Springer, New York, NY. https://doi.org/10.1007/978-1-4419-0118-7.

Robinson MD, Smyth GK (2008). Small sample estimation of negative binomial dispersion, with applications to SAGE data. Biostatistics 9, 321-332. https://doi.org/10.1093/biostatistics/kxm030

Rainey C (2024). Hypothesis tests under separation. Political Analysis 32, 172-185. https://doi.org/10.1017/pan.2023.28

Greenland S (1992). Likelihood-ratio testing as a diagnostic method for small-sample regressions. Annals of Epidemiology 2, 311-316.