is this odds ratio transformation correct?

Question

I am using a function based on the r package stargazer to transform outcomes from a poisson analysis to odds ratio. I want to make sure that the output is correct, because I don't fully trust the function I used.

Here are the summary data from the regression models:

summary(TNC_count)
Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-51.94  -18.21  -13.93   -2.51  121.78  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) 5.107e+00  3.077e-03  1659.7   <2e-16 ***
NFIP_count  4.401e-03  1.563e-05   281.6   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 335093  on 633  degrees of freedom
Residual deviance: 294111  on 632  degrees of freedom
AIC: 296801

Number of Fisher Scoring iterations: 6



summary(ncdem_model)
Call:
glm(formula = NCDEM_count ~ NFIP_count, family = "poisson", data = dmg_model)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-95.724  -21.345  -17.824    0.185  100.932  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) 5.423e+00  2.599e-03  2086.6   <2e-16 ***
NFIP_count  5.068e-03  1.083e-05   468.1   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 466864  on 633  degrees of freedom
Residual deviance: 361648  on 632  degrees of freedom
AIC: 364279

Number of Fisher Scoring iterations: 6

and here is how it looks after using stargazer to output it as odds ratio, with (correct me if I'm wrong) the coefficients in odds ratio with the standard errors underneath.

stargazer2(list(tnc_model, ncdem_model), odd.ratio = T, type = "text")

==============================================
                      Dependent variable:     
                  ----------------------------
                    TNC_count     NCDEM_count 
                       (1)            (2)     
----------------------------------------------
NFIP_count           1.004***      1.005***   
                    (0.00002)      (0.00001)  
                                              
Constant            165.156***    226.647***  
                     (0.508)        (0.589)   
                                              
----------------------------------------------
Observations           634            634     
Log Likelihood     -148,398.700  -182,137.700 
Akaike Inf. Crit.  296,801.400    364,279.500 
==============================================
Note:              *p<0.1; **p<0.05; ***p<0.01
> stargazer2(list(nfip_model1, nfip_model2, nfip_model3), odd.ratio = T, type = "text")

Is this correct?

For what it's worth, this is the function I am using with stargazer:

stargazer2 <- function(model, odd.ratio = F, ...) {
  if(!("list" %in% class(model))) model <- list(model)
  
  if (odd.ratio) {
    coefOR2 <- lapply(model, function(x) exp(coef(x)))
    seOR2 <- lapply(model, function(x) exp(coef(x)) * summary(x)$coef[, 2])
    p2 <- lapply(model, function(x) summary(x)$coefficients[, 4])
    stargazer(model, coef = coefOR2, se = seOR2, p = p2, ...)
    
  } else {
    stargazer(model, ...)
  }
}

Answer 1

This is incorrect, primarily because exponentiating the coefficients in a Poisson model does not give you an odds ratio. What it would be the odds of? Exponentiating the coefficient gives you an incidence rate ratio (IRR). You can also use Poisson models with binary outcomes, in which case the exponentiated coefficients are interpreted as risk ratios.

The second thing I see incorrect is the formula for computing the standard error of the IRR. This value should not be computed; it does not provide any useful information because it cannot be used to describe symmetric variation in the IRR or be used in a z-test of the IRR. The standard error of the log IRR (i.e., of the coefficients) can be used in these ways and would be appropriate to report along with the (un-exponentiated) coefficients. You can, however, report confidence interval bounds on the IRR scale by exponentiating the coefficient confidence interval bounds, and these might make more sense to report.