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, ...)
}
}
```

## Best Answer

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.