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.

Can anyone help/advise ?

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.

  • $\begingroup$ Does removing the intercept from the model with - 1 help? glm(cases ~ agecat - 1, family = poisson(link = "log") $\endgroup$ – jthetzel Nov 16 '11 at 19:06
  • $\begingroup$ 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. $\endgroup$ – P Sellaz Nov 16 '11 at 19:22
  • $\begingroup$ 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. $\endgroup$ – jthetzel Nov 16 '11 at 19:34
  • $\begingroup$ 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. $\endgroup$ – P Sellaz Nov 16 '11 at 19:38
  • $\begingroup$ Great. See: tinyurl.com/reproducible-000 $\endgroup$ – jthetzel Nov 16 '11 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).

| cite | improve this answer | |
  • $\begingroup$ I tried this already:glm(formula = cases ~ agecat,offset=log(population), family = poisson(link = "log")) but it doesn't give the correct result. The effect of exp(pop) on the stata model is very small, but the effect of offset=log(population) in the R model is huge (after exponentiation) $\endgroup$ – P Sellaz Nov 16 '11 at 19:16
  • $\begingroup$ 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. $\endgroup$ – onestop Nov 16 '11 at 19:56

Have you looked at glm documentation? there is an offset() function that can be included.

glm(cases~agecat+offset(log(population)), family = poisson(link = "log"))
| cite | improve this answer | |
  • $\begingroup$ Yes, I tried offset=log(population) but it doesn't give the correct result. $\endgroup$ – P Sellaz Nov 16 '11 at 19:20
  • $\begingroup$ @P Sellaz: It appears that @B_Miner's suggestion is for a function that goes in the glm formula, as opposed to an argument to glm. I'd guess they act the same, but perhaps not. $\endgroup$ – Wayne Nov 16 '11 at 19:29
  • $\begingroup$ They seem to act the same $\endgroup$ – P Sellaz Nov 16 '11 at 19:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.