Calculate interaction effect confidence intervals in zero-inflated poisson regression

I'm conducting a zero-inflated Poisson regression using the pscl package in R. I've included interaction terms but am having an issue with interpretation. I am assuming an additive effect and summing coefficients (x + y + xy) but am not sure what to do about the confidence intervals or p-values. From what I understand, these have to be re-estimated, probably with some bootstrapping method, but I can't figure out how to do this.

The main issue is that one of my effects reverses in interaction. I've provided a simplified version of the code below (sorry, I can't share the data). Here's a brief description of the scenario: when doctors discuss re-injury prevention with their patients, time off work increases, but in interaction with a low-stress setting, it reduces time off work.

So my question is:

• How do you calculate interaction effect confidence intervals and p-values using a zeroinfl object?
• Is the process for calculating CIs and p-values different between the zero and count models?

Example code would be greatly appreciated!

model <- zeroinfl(TimeLoss ~ PrevDisc + LowStress +
PrevDisc * LowStress,
data = Doctors)


If PrevDisc and LowStress are binary variables (this is my impression from your description), then the interaction model simply corresponds to four different zero-inflated Poisson distributions: one for each combination of PrevDisc and LowStress.

When using the formula TimeLoss ~ PrevDisc * LowStress, treatment coding of the coefficients is used, i.e., the four parameters (in each model part) are coded as an intercept, two main effects, and an interaction effect. This coding facilitates judging whether or not the interaction effect is significant.

If you want to assess the PrevDisc effect separately for the two LowStress groups, then you can use a nested coding of the coefficients via the formula TimeLoss ~ LowStress/PrevDisc.

In either case, all inference can be done "in the usual way", i.e., summary() and confint() for marginal Wald tests and Wald confidence intervals. But also lrtest() (from lmtest) for nested model comparisons with the likelihood ratio test or AIC()/BIC(). (Generalized) linear hypotheses can be tests with linearHypothesis() (from car) and glht() (from multcomp) respectively.

The difference in interpretation between the two model parts is that the count model is a log-linear model for the mean in the count component. The zero model is log-linear for the odds of zero inflation (i.e., an observation from the point mass component).

• Thanks for the response but I'm not sure I understand. Does this mean the 'summary()' and 'confint()' functions will give me the figures I need and there's no need to bootstrap? Assuming I'm looking at the raw coefficients, I could sum the main effects and interaction effect to get the total effect; does that mean I can use 'confint()' to get the interaction CI and just adjust it based on the difference between interaction effect and total effect? I.e., if +0.5, just increase the CIs by +0.5 each? And then from there see whether the confidence intervals cross the significance threshold? – Tyler Lane Feb 20 '19 at 5:52
• What exactly do you want to have confidence intervals for? – Achim Zeileis Feb 20 '19 at 6:28
• The effect reverses in interaction: one of the main effects is negative, but the interaction effect positive. What I am trying to figure out is whether I can interpret this reversal as neutralising the negative main effect, or it actually becomes positive. CIs for the fully summed interaction effect (x + y + xy) should tell me this, I think. However, I am open to the possibility that I've completely misunderstood my problem. – Tyler Lane Feb 20 '19 at 20:23
• I'm still not sure what exactly you are looking for. Do you want a nested coding of the PrevDisc effect, separately for each LowStress group? Then you can use the formula TimeLoss ~ LowStress/PrevDisc instead of the usual treatment effect coding with *. I have modified my answer to explicitly incorporate this. – Achim Zeileis Feb 24 '19 at 1:04
• Thanks for your assistance. I don't think that's what I'm looking for. There's a function in STATA that achieves this, 'lincom', which "computes point estimates, standard errors...p-values, and confidence intervals for linear combinations of coefficients after any estimation command." stata.com/manuals13/rlincom.pdf. I'm speaking with a biostatistician in my department this afternoon and will update this post if we arrive at a solution. – Tyler Lane Feb 24 '19 at 22:43

I found the biostat3 that applies STATA's lincom function, which works with zeroinfl models. So, given the zeroinfl model:

model <- zeroinfl(TimeLoss ~ PrevDisc + LowStress +
PrevDisc * LowStress,
data = Doctors)


You can use the following function to derive point estimates, confidence intervals, and p-values of the interaction effect. I was interested in an interaction effect specifically affecting the count model so report that one below, hence the "count_" prefix; the zero-models can be derived using "zero_" prefix.

lincom(model, c("count_PrevDisc + count_LowStress + count_PrevDisc:LowStress",
form = TRUE)].


The responses are exponentiated, meaning you get the Odds Ratios for the zero model and Incidence Relative Risk right away.