1

Im trying to plot the marginal effect of a specific variable in a poisson regression and then correlate that graphic with its corresponding incidence rate ratio.

I've achieved this for most of my plots. However for one of them, the incidence rate ratio signals an overall positive association of my variable of interest and the plot shows a clearly negative association. From my understanding, there should be something wrong with this.

Could you help me ? :) I might be understanding something wrong in my analysis...

I am first creating the Poisson model:

model3<- glm(y ~ x1*x2 + x3 + x4 + x5, data=data, family = poisson)

From which I get the following IRRs

poissonirr(y ~ x1*x2 + x3 + x4 + x5, data=data)

Incidence-Rate Ratio:
                         IRR  Std. Err.        z     P>|z|    
x1                 1.03404133 0.00471847   7.3359 2.202e-13 ***
x2                 1.16795382 0.01235611  14.6752 < 2.2e-16 ***
x3                 0.63214010 0.00817795 -35.4523 < 2.2e-16 ***
x4                 1.00468920 0.00095329   4.9305 8.204e-07 ***
x5                 0.98118299 0.00267124  -6.9776 3.003e-12 ***
x1:x2              0.99382845 0.00073716  -8.3462 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Then I plot the marginal effect of the first variable in the model (x1) and I get the following plot:

plot_model(model3, type = "eff", terms = c("x1"))

x1 vs y Poisson Regression plot (with negative association)

which clearly shows an apparent negative association between x1 and y

thank you in advance !!

(I am using the mfx package to calculate IRRs and sjPlot::plot_model for the plotting)

Miquel Serna Pascual
  • 11
  • 2
  • You might get help from someone very familiar with this package (you don't mention - is it mfx?) or with incident rate ratios. But without example data or the code used to produce the plot, not many can answer. – markhogue Nov 08 '19 at 13:26
  • Thanks @markhogue, I just posted the code I used for the plotting and yes, the package I am using for IRR is mfx and for plotting plot_model – Miquel Serna Pascual Nov 08 '19 at 13:32
  • So, just an observation, and I don't know how anything about this - could the problem be that the model is x1*x2, not x1 by itself? x1*x2 has a negative slope, but both x1 and x2 by themselves show positive. Could be plot_model is not grabbing the model details you expect. – markhogue Nov 08 '19 at 14:05
  • It might have something to do with the interaction, is it possible to share a subset of your data? Otherwise it's quite hard to figure out what's so weird about this regression. – StupidWolf Nov 08 '19 at 23:36
  • And just to add on, the incidence rate is just the exponential of the coefficients from poisson regression, hence irr < 1 means a negative coefficient under poisson regression – StupidWolf Nov 08 '19 at 23:37
  • Thank you very much everyone! Indeed, it was an issue with the interaction. The IRR refers to the effect of x1 when x2 is 0, and the effect changes as described by the IRR of the interaction itself. Thus, the IRR was not representative of the effect of x1 in the whole dataset. The plot was right. The solution we found was to plot the effect of x1 at different values of x2 (via the "int" option of plot_model, so we got a much better picture of what was going on and we successfully reached our conclusions. Thanks @markhogue and @StupidWolf! Really appreciate it! :) – Miquel Serna Pascual Nov 12 '19 at 08:54

1 Answers1

0

Since an interaction is involved, you can't actually just interprete the main effect alone, but need to take the effect of the interaction into account. Thus, I would recommend plotting the following:

plot_model(model3, type = "eff", terms = c("x1", "x2"))

The effect of x1 is different for different values / levels of x2. Thus, your plot that only considers x1 is misleading (and so is the coefficent of x1, when you don't also look at the coefficient for x1:x2).

Daniel
  • 7,252
  • 6
  • 26
  • 38