6
$\begingroup$

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.

$\endgroup$
6
  • $\begingroup$ Hi newspice! What are you modeling exactly? Number of seeds germinated? Probability of seeds germinated? How many seeds per treament? Can share a subset of your data structure? What does this exactly mean "For example the below is the summary output of the model when only filtered to one seed with around 24 observations." $\endgroup$
    – Stefan
    Mar 21 at 23:07
  • $\begingroup$ Thanks Stefan. Ideally, I would like to model the probability that seeds will germinate but really I'm just looking to model the effect of treatments on germination so the response variable is flexible. I have edited my question to include the sample data I used in the example model outputs. My overall dataset has 11 different seed species so to see the effect of treatments on an individual seed I filter the dataset to just 1 seed. $\endgroup$
    – newspice
    Mar 21 at 23:30
  • $\begingroup$ Thanks! When you fit the glmm, what was your random effect? Why did you move away from it? Also why wouldn't you want to add species into the model as a fixed effect? The problem you are describing is called complete separation, i.e. when you have all zeros or ones for a given treatment combination. See here and here. There are many more examples here on CV regarding this. $\endgroup$
    – Stefan
    Mar 21 at 23:51
  • 1
    $\begingroup$ What do you think about analyzing this as a two step model, i.e. a hurdle model? First you ask what is the probability of germination and then using a zero-truncated poisson, what's rate of germination for those seeds that germinated? All of this can be done in one model using the glmmTMB package. $\endgroup$
    – Stefan
    Mar 21 at 23:57
  • 1
    $\begingroup$ Thanks so much Stefan. I wasn't aware there was a term for it! My random effect was originally the species for the overall model fitted by adding (1 | species) however for an individual species set like the example, I thought this was unnecessary to include given the species remains the same for the entire set. I will look to implementing a two step model like you mentioned and will post my results. $\endgroup$
    – newspice
    Mar 22 at 0:40

1 Answer 1

8
$\begingroup$

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.

$\endgroup$
1
  • 2
    $\begingroup$ Thanks so much Gordon. That makes a lot of sense and I was able to implement this across other seeds in my dataset successfully. Much appreciated! $\endgroup$
    – newspice
    Mar 22 at 11:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.