* Stata
poisson y x1 x2, irr vce(robust)          // Poisson + 稳健标准误
nbreg y x1 x2, irr                        // 负二项回归
poisson y x1 x2, irr exposure(t)          // 随访时间不等：按 t 调整

## R
glm(y ~ x1 + x2, data = df, family = poisson)
library(MASS)
fit <- glm.nb(y ~ x1 + x2, data = df)     // 负二项
exp(cbind(IRR = coef(fit), confint(fit)))