It seems that relative risk estimation using Poisson regression is a popular approach in epidemiology, with a couple of papers I have recently looked at (for example [3] and many others like it) following an approach outlined by G. Zou in [1]. However, I seem to have some difficulties following the logic of the approach and probably misunderstood something, so I hope someone can shed some light on this for me.
For binary outcome y (e.g. contracting a disease) and binary exposure x, Zou proposes we set
$\log \pi(x_i) = \alpha + \beta x_i$
and refers to $\pi(x_i)$ as "an underlying risk" for subject $i$, which seems to suggest that it should be the probability of $y_i=1$, i.e., $Pr[y_i = 1|x_i] = \pi(x_i)$, based on common interpretation of the term "risk" (as explained in eg. [2]). This view is supported by the claim in [1] that then the relative risk is
$RR = \frac{Pr[y_i = 1|x_i = 1]}{Pr[y_i = 1|x_i = 0]} = \frac{\pi(x_i=1)}{\pi(x_i=0)} = \frac{e^{\alpha+\beta}}{e^\alpha}= e^\beta .$
However, Zou then fits this using a Poisson log-likelihood using $\pi(x_i)$ as the Poisson rate parameter $\lambda$, i.e., $C+\sum_i [y_i \log(\pi(x_i) - \pi(x_i)]$
which of course assumes that $Pr[y_i|x_i]$ follows a Poisson distribution and is therefore not the same as $\pi(x_i)$.
Now, assuming the $Pr[y_i|x_i]$ is Poisson, we get as relative risk
$RR = \exp( \log( \frac{Pr[y_i = 1|x_i = 1]}{Pr[y_i = 1|x_i = 0]} )) = \exp( \alpha + \beta - e^{\alpha+\beta} - (\alpha - \exp(\alpha)) ) = e^\beta \exp(e^{\alpha+\beta}-e^\alpha)$
Similar calculations would lead us to conclude that $\exp(\beta)$ is actually the value for the odds ratio:
$odds(x_i) = \exp( \log( \frac{Pr[y_i = 1|x_i]}{Pr[y_i = 0|x_i]} )) = \exp( \alpha+\beta x_i - e^{\alpha+\beta x_i} + e^{\alpha+\beta x_i} ) = e^{\alpha+\beta x_i} = \pi(x_i),$
$OR = \frac{odds(x_i=1)}{odds(x_i=0)} = \frac{\pi(x_i=1)}{\pi(x_i=0)} = e^\beta.$
It seems to me like something does not add up here - Zou's Poisson regression approach is used in the papers referenced above - which explicitly use this to obtain the relative risk and NOT an odds ratio - and implemented in popular statistical software packages (such as Stata and R), none of which acknowledges the above apparent contradiction. Given that, it seems likely that my logic is flawed, but I cannot see where and would appreciate any help.
There is this answer to a question about Poisson regression and OR which states that "exponentiating the coefficients gives you the multiplier corresponding to a unit change in the predictor variable. [...] It is not an odds ratio." which is in contradiction to my calculations above, but that might be due to the special case of binary x_i here whereas the quoted answer above had no such restriction.
There is also the question Poisson regression to estimate relative risk for binary outcomes which asks why the approach is not used more often. In the accepted answer it is pointed out that " A Poisson regression is estimating often a rate, not a risk, and thus the effect estimate from it will often be noted as a rate ratio (mainly, in my mind, so you can still abbreviate it RR) or an incidence density ratio (IRR or IDR)" - so is this all just down to incorrect use of the term "risk"?
[1]: G. Zou: A Modified Poisson Regression Approach to Prospective Studies with Binary Data, https://academic.oup.com/aje/article/159/7/702/71883
[2]: C. Andrade: Understanding Relative Risk, Odds Ratio, and Related Terms: As Simple as It Can Get, https://www.pitt.edu/~bertsch/risk.pdf
[3]: C. Niedzwiedz et al.: Ethnic and socioeconomic differences in SARS-CoV-2 infection: prospective cohort study using UK Biobank, https://bmcmedicine.biomedcentral.com/articles/10.1186/s12916-020-01640-8