Age-adjusted rates using a poisson regression

[Temi Schultz]

Hello all. I'm using a poisson regression calculate age-adjusted rates and incidence rate ratios for a large dataset. I'd like to calculate age-adjusted rate for specific sub-populations. Eg. the age-adjusted rate incidence rate for current smokers. However, I do not know how to do this without stata also adjusting for smoking status. The codes I'm using at the moment are:

poisson _d i.agegrp i.smoke, exp(pyear) irr

margins, at(smoke==3) predict(ir)

I believe that in typing i.smoke, stata automatically also adjusts for smoking status. However if I do not include smoking in the regression, then stata is not able to calculate the margins for that group. Is there a way to get the value for current smokers, without adjusting for smoking status in the process of doing so?


Thank you in advance for your help!

[Temi Schultz]

Thanks for the addition Andrea, I have now made this edit.

Best,
Temi

[Andrea Discacciati]

To follow up on Clyde's excellent answer.

If what you're trying to do is to reproduce with -poisson- the results you obtained from non-parametric direct standardisation (-dstdize-), then you need to include interaction terms between agegrp and smoke in your Poisson model.

See the following example, where I directly sex-standardise the incidence rate of the outcome by levels of stage. The Poisson regression without the interaction makes modelling assumptions that the non-parametric standardisation doesn't make, hence the difference in standardised IRs.

. dstdize status01 surv_mm sex, by(stage) using(`w')

[...output omitted...]

Summary of Study Populations:
    stage             N      Crude     Adj_Rate       Confidence Interval
 --------------------------------------------------------------------------
        0        124021   0.004507     0.004555    [  0.004177,    0.004933]
        1        466178   0.003863     0.003941    [  0.003758,    0.004124]
        2         18178   0.014303     0.014223    [  0.012436,    0.016009]
        3         10747   0.040476     0.040976    [  0.037034,    0.044918]

.
. set cformat %7.6f

. qui poisson status01 i.stage i.sex, exposure(surv_mm)

. margins stage,predict(ir)

Predictive margins                              Number of obs     =      7,775
Model VCE    : OIM

Expression   : Predicted incidence rate, predict(ir)

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       stage |
    Unknown  |   0.004557   0.000193    23.63   0.000     0.004179    0.004935
  Localised  |   0.003918   0.000093    42.27   0.000     0.003736    0.004100
   Regional  |   0.013824   0.000860    16.07   0.000     0.012138    0.015510
    Distant  |   0.039044   0.001883    20.73   0.000     0.035353    0.042735
------------------------------------------------------------------------------

.
. qui poisson status01 i.stage##i.sex, exposure(surv_mm)

. margins stage, predict(ir)

Predictive margins                              Number of obs     =      7,775
Model VCE    : OIM

Expression   : Predicted incidence rate, predict(ir)

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       stage |
    Unknown  |   0.004555   0.000193    23.55   0.000     0.004176    0.004934
  Localised  |   0.003941   0.000093    42.17   0.000     0.003758    0.004124
   Regional  |   0.014223   0.000918    15.49   0.000     0.012424    0.016022
    Distant  |   0.040976   0.002054    19.95   0.000     0.036950    0.045002
------------------------------------------------------------------------------

[Temi Schultz]

Thank you very much for the clarification, I followed this and have now achieved my desired output! This was really clear and hugely useful.

Best wishes,
Temi

[Clyde Schechter]

No, it is not smoking-adjusted. Here's how it's computed (more or less--the computational algorithm inside -margins- may be different, but this is what it computes):

1. Start with the base level of the smoke variable. Let's say that's 1 just for concreteness. Set smoke = 1 in the entire data set.

2. Set smoke equal to 1 in the entire data set.

3. Calculate the predicted incidence rate in each observation of the data set.

4. Average the results from step 3 and report that as the average incidence rate for smoke = 1.

5. Repeat steps 2 through 4, using smoke = 2, and smoke = 3, ...

As you can see, this calculation does adjust for age, because the ages in the data are used as is in each calculation. But it does not adjust for smoking because in each cycle, only one value of smoke is used.

This comports with the standard definition of average age-adjusted incidence rate, and it is what I (and any other epidemiologist) would expect to see for a statistic reported under that name.

[Temi Schultz]

Thank you Clyde. I have done this also, but I believe that the output of the code you suggested is then age AND smoking-adjusted rate, as opposed to purely age-adjusted rate. My thought is that the above code you suggested is the same method adopted when adjusting for several variables. However do correct me if I am wrong.

Best,
Temi

[Clyde Schechter]

If I understand correctly what you want to do, it would be:

poisson _d i.agegrp i.smoke, exp(pyear) irr
margins smoke,predict(ir)

This will do a full-sample adjustment for age and give you the age-adjusted incidence rates in all categories of the smoke variable.

[Temi Schultz]

Thank you for your swift response Carlo. I have tried this but the problem is that in doing this Stata only adjusts age for the subpopulation of current smokers.
However, I would like age adjusted across the whole population in order to then compare the age-adjusted rates for current smokers with that of past smokers and never smokers.

Best,
Temi

[Carlo Lazzaro]

Waiting for your username to be fixed, you may want to try:

poisson _d i.agegrp if smoke==3, exp(pyear) irr

[Temi Schultz]

Apologies, I am a new member and overlooked this upon registration. I have now contacted administration to rectify it.

Many thanks for letting me know.

[Carlo Lazzaro]

As per FAQ, please note that real given and family names are preferred on this forum. Thanks.