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!

[Carlo Lazzaro]

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

[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]

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

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

[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

[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 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]

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 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

[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]

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

Best,
Temi

[Md Bayzidur Rahman]

Can someone please give me a reference of the adjustment method Clyde described in 5 steps?

[Clyde Schechter]

I don't think you can ever get a reference to support a particular snippet of code. But if you want a reference about the -margins- command and how it does adjusted predictions, I recommend: Williams R. Using the margins command to estimate and interpret adjusted predictions and marginal effects. The Stata Journal (2012) 12(2) 308-331.

[Md Bayzidur Rahman]

Hi Clyde,
Thanks for the reference. I was looking for a stats reference for calculating average covariate adjusted rates. I can see that different people does this in different ways. For example, I found some one estimated the adjusted rates in the following way
Step 1: Fitted a GEE poisson model for the panel variable including the potential adjustment variables in the model
Step 2: predicted the counts for all the observations in the dataset
Step 3: Aggregated the predicted counts and the observed counts over panels
Step 4: calculated the overall observed rate as Total observed counts in all the panels/ Total days at risk
Step 5: calculated the adjusted panel specific rate = panel specific observed counts (from step 3)/ panel specific expected counts (from step 3)*overall observed rate ( Step 4)
Step 6: Estimated the CI of the adjusted rates by doing manual bootstrapping
I couldn't find any methodological reference for this approach. I also prefer the margins approach as you explained, which is similar to directly adjusted survival curves from a PH model.

[Clyde Schechter]

Well, as you have pointed out, there are many ways to adjust predicted values. None of them can be distinguished as "the right way" to do it. I think the best one can hope for is an explanation of the methodology used, in transparent enough terms, that you can figure out what was done and what it means and what conclusions and comparisons can be reasonably based on it. The -margins- command is widely used in the Stata community. And other major statistical packages also have commands that do the same.

One reference which proposes this method for use in clinical studies is Chang I-M, Gelman R, Pagano M. Corrected group prognostic curves and summary statistics. J Chron Dis 1982 35:669-674. I do not know if the method originates with these authors, or if they are merely introducing it into the medical literature.