CIs in predictions from zero-inflated negative binomial model

[martinspn]

Hi all

I am fitting a zero-inflated model to estimate a count that ranges from 0 to 10 from a similar variable that also ranges between these values. I may end up using an ordinal or a zero-inflated beta regression model due to some values being in the upper limit. However, I came across an intriguing result when fitting the negative binomial model and would like to hear your insights.

I have two questions:

1. When estimating the probability of outcome == 0, why is lower bound of the CI for the estimate when aspects_score_man == 0 negative?
2. It seems to me that the CI bounds are being calculated as the estimate ± 1.96*SE using the standard errors in the response scale. Shouldn't these be calculated in the link scale and then exponentiated? I usually do that for logistic regression models, but I did not find such possibility for this case.

Below is the code I used to fit the model:

summary(m1 <- zeroinfl(aspects_staining_rev ~ rcs(aspects_score_man, c(5, 7, 8)),
                                         data = df_xper_final, dist = "negbin"))

Call:
zeroinfl(formula = aspects_staining_rev ~ rcs(aspects_score_man, c(5, 7, 8)), data = df_xper_final, dist = "negbin")

Pearson residuals:
    Min      1Q  Median      3Q     Max 
-1.1435 -0.6774 -0.5203  0.4621  6.6918 

Count model coefficients (negbin with log link):
                                                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)                                           2.27279    0.27305   8.324  < 2e-16 ***
rcs(aspects_score_man, c(5, 7, 8))aspects_score_man  -0.20252    0.04779  -4.238 2.26e-05 ***
rcs(aspects_score_man, c(5, 7, 8))aspects_score_man'  0.06312    0.03863   1.634    0.102    
Log(theta)                                            1.50028    0.26696   5.620 1.91e-08 ***

Zero-inflation model coefficients (binomial with logit link):
                                                     Estimate Std. Error z value Pr(>|z|)  
(Intercept)                                          -1.29572    0.75755  -1.710   0.0872 .
rcs(aspects_score_man, c(5, 7, 8))aspects_score_man   0.11787    0.12593   0.936   0.3493  
rcs(aspects_score_man, c(5, 7, 8))aspects_score_man'  0.09753    0.08789   1.110   0.2671  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Theta = 4.4829 
Number of iterations in BFGS optimization: 36 
Log-likelihood: -1524 on 7 Df

And below are the estimated probabilities of zero for each value of the predictor.

> marginaleffects::predictions(m1,
+                              newdata = data.frame(aspects_score_man = 0:10),
+                              type = "zero")

 Estimate Std. Error     z Pr(>|z|)     S   2.5 % 97.5 % aspects_score_man
    0.215     0.1278  1.68  0.09270   3.4 -0.0356  0.465                 0
    0.235     0.1142  2.06  0.03921   4.7  0.0116  0.459                 1
    0.257     0.0979  2.63  0.00859   6.9  0.0654  0.449                 2
    0.280     0.0793  3.54  < 0.001  11.3  0.1250  0.436                 3
    0.305     0.0592  5.15  < 0.001  21.9  0.1888  0.421                 4
    0.330     0.0404  8.17  < 0.001  51.5  0.2512  0.410                 5
    0.359     0.0322 11.17  < 0.001  93.8  0.2964  0.423                 6
    0.405     0.0332 12.21  < 0.001 111.4  0.3401  0.470                 7
    0.477     0.0259 18.39  < 0.001 248.5  0.4260  0.528                 8
    0.555     0.0248 22.37  < 0.001 365.8  0.5062  0.603                 9
    0.630     0.0336 18.78  < 0.001 258.9  0.5645  0.696                10

Columns: rowid, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, aspects_score_man, aspects_staining_rev 
Type:  zero 

CI calculations by hand:

> 0.215 + 1.96*0.1278
[1] 0.465488
> 0.215 - 1.96*0.1278
[1] -0.035488

These values match the ones obtained from predictions().

Then I discussed with @arthur-albuquerque and decided to use bootstrapping to see what would change. The CIs make more sense, having bounds that are positive and asymmetric relative to the estimate in the response scale.

> marginaleffects::predictions(m1,
+                              newdata = data.frame(aspects_score_man = 0:10),
+                              type = "zero") %>% inferences(method = "boot")

 Estimate Std. Error  2.5 % 97.5 % aspects_score_man
    0.215     0.1330 0.0400  0.542                 0
    0.235     0.1185 0.0586  0.516                 1
    0.257     0.1021 0.0865  0.490                 2
    0.280     0.0837 0.1240  0.464                 3
    0.305     0.0634 0.1786  0.431                 4
    0.330     0.0431 0.2368  0.414                 5
    0.359     0.0319 0.2950  0.421                 6
    0.405     0.0320 0.3421  0.468                 7
    0.477     0.0250 0.4271  0.527                 8
    0.555     0.0244 0.5076  0.602                 9
    0.630     0.0335 0.5589  0.692                10

Columns: rowid, estimate, aspects_score_man, aspects_staining_rev, std.error, conf.low, conf.high 
Type:  zero 

> 0.215 + 1.96*0.1330
[1] 0.47568
> 0.215 - 1.96*0.1330
[1] -0.04568

I'd love to know your thoughts on why this is happening and if I did anything inappropriate when getting the estimates.

Thanks!!

[vincentarelbundock]

Yes, your guess is right: intervals are built symmetrically around the estimates.

We can only back transform in a few models like glm() when the predict() method allows to making predictions on the link scale and the inverse function is known.

The pscl package doesn't seem to have a "zerolink" option (or similar), so I don't think this is possible

https://www.rdocumentation.org/packages/pscl/versions/1.5.9/topics/predict.zeroinfl

[martinspn]

Thank you @vincentarelbundock!