# Incidence rate ratios (Stata vs R)

I asked this question on Stackoverflow: https://stackoverflow.com/questions/8142118/incidence-rate-ratios-in-r-poisson-regression and was advised to post here instead. I have data that looks like this

   sex agecat  cases population

1 male    0-4  12     126526
2 male    5-9  12     128375
3 female  0-4  11     129280
4 male    10-14 4     127910
5 female  0-4  13     127158
6 male    0-4  8      125125


I want to duplicate the output I get in stata with this command

poisson cases i.agecat, exp(pop) irr


which gives output such as:

   cases |        IRR   Std. Err.      z    P>|z|     [95% Conf. Interval]
---------+----------------------------------------------------------------
agecat |
2  |   .5125755   .0530442    -6.18   0.000     .4578054    .6669639
3  |    .323456   .0381304    -9.60   0.000     .2665044    .4172274
population | (exposure)


in R with a command such as

glm(cases~agecat, family = poisson(link = "log")


I know I need to exponentiate the coefficients and confidence intervals, but I think I also need some kind of offset so so that the intercept is zero; and adjust for per unit population vs baseline.

Thanks EDIT: The question on SO has been answered, but I posted more detail here. In particular, I think the issue has to do with adjusting for population size in stata with exp(pop) - and how to replicate this in R.

• Does removing the intercept from the model with - 1 help? glm(cases ~ agecat - 1, family = poisson(link = "log") Nov 16, 2011 at 19:06
• No, this just puts the estimate for the intercept to the estimate of the baseline effect of agecat. I think the intercept issue is a red herring actually. It seems that using the irr option in stata just suppresses the display of the intercept. Nov 16, 2011 at 19:22
• As Ben mentioned in SO, it would be much easier to help if you provided a reproducible example, along with your result from Stata and R. Nov 16, 2011 at 19:34
• Hi jthetzel, yes, understood. I'm making a small reproducible dataset. What is the best way to post it ? It's just 9 rows similar to the above data. Nov 16, 2011 at 19:38
• Great. See: tinyurl.com/reproducible-000 Nov 16, 2011 at 19:45

The Stata option exp(pop) includes log(pop) as an offset term in the linear predictor, so the R equivalent should be offset=log(pop).
• Adding offset=log(population) to the R model should make a big difference to the estimate of the constant, but little difference to the estimates of the other parameters, as the populations are all similar. If it makes a big difference to those as well, something's going wrong. Nov 16, 2011 at 19:56
glm(cases~agecat+offset(log(population)), family = poisson(link = "log"))