# Non-intuitive answer from a Poisson regression

I am not very familiar with Poisson regression, so I think I must have made a mistake in the below analysis:

I am studying the effects of smoking on lung cancer rates. The dataset is provided here. The variable smoking_status is defined:

smoking status: coded 1 = doesn't smoke, 2 = smokes cigars or pipe only, 3 = smokes cigarrettes and cigar or pipe, and 4 = smokes cigarrettes only,

I modified the data a bit and made two new categorical variables: pipe/cigar smoker and cigarette smoker, to replace smoking_status. So a smoking status 1 maps to (0,0), 2 maps to (1,0), 3 maps to (1,1) etc.

I also added a constant column to my dataset. This is all I did to the data.

I then performed Poisson regression on this dataset, using a exponential link function. My hopes was that the coefficients of the two new variables would be positive, but instead only cigarette_smoker is positive. The confidence intervals do not contain positive points either.

Have I analysed the data incorrectly, or is my data just wrong?

EDIT

The output (it's from a Python library Statsmodels )

             Generalized Linear Model Regression Results
Dep. Variable:                      y   No. Observations:                   36
Model:                            GLM   Df Residuals:                       31
Model Family:                 Poisson   Df Model:                            4
Link Function:                    log   Scale:                             1.0
Method:                          IRLS   Log-Likelihood:                -815.93
Date:                Thu, 31 Jan 2013   Deviance:                       1391.8
Time:                        13:19:32   Pearson chi2:                 1.22e+03
No. Iterations:                     7

coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
x1             0.2596      0.006     44.097      0.000         0.248     0.271
x2            -0.1850      0.024     -7.775      0.000        -0.232    -0.138
x3             0.5327      0.031     17.101      0.000         0.472     0.594
x4             0.0004   7.95e-06     54.637      0.000         0.000     0.000
const          2.9593      0.046     63.903      0.000         2.869     3.050


The variables in order are age, smoke_cigar (0,1), smoke_cigarettes (0,1), population (in hundred of thousands), constant_term.

Some example data:

array([[ 2., 0., 0., 359., 1.], [ 4., 0., 1., 3270., 1.]])

with target deaths [ 22., 514.] respectively.

• Pipe/cigar smokers don't typically inhale the smoke, or not much of it. You're only looking at lung cancer, so perhaps it's not so counter-intuitive. – Scortchi - Reinstate Monica Jan 31 '13 at 15:53
• mmmm but it produces a negative coefficient; surely it can't prevent lung cancer. – Cam.Davidson.Pilon Jan 31 '13 at 16:11
• Hi, Cam. It would be helpful to see some output to potentially help diagnose things. Note that what you've done is recoded things by creating two main effects and effectively removing an interaction (though the model is equivalent). Also, are you controlling for exposure or just modeling raw counts? If the latter, that would probably explain it: There just aren't many pipe/cigar smokers compared to nonsmokers and cigarette smokers in the first place. – cardinal Jan 31 '13 at 17:46
• My first thought (before fully reading cardinal's note above, which covered my points) was that the original model was treating cigar-smoking and cigarette-smoking as independent (rather than including an interaction term)? Also, I'm not sure about how one specifies an offset (i.e. denominator) variable in Poisson regression in Python: usually this wouldn't show up with its own coefficient in any of the stats packages I use (SAS, Stata, R) -- see StasK's answer to this question for info on the offset variable stats.stackexchange.com/questions/23672/… – James Stanley Jan 31 '13 at 20:36
• I think the main issue here turns out to be non-specification of the offset variable (which @cardinal spotted this) so that you were modelling counts rather than rates. I've discussed this and the interaction issue in the answer below: it's rather longer than intended, mostly because I'm currently revising materials for a 2-day course on regression methods in a couple of weeks and adding in lots of R code to that (so this felt like a legitimate use of my time...) – James Stanley Jan 31 '13 at 23:06

## 1 Answer

A brief (actually it turned out quite lengthy) update on material in the comments (as covered by @cardinal and myself above):

I think the main problems are arising from not specifying an offset variable (as per my comment above, see StasK's note on Interpretation of Coefficient answer) which would usually be the natural log of the population at risk (or person-time at risk.)

For Poisson regression to compare rates rather than counts, you need to specify an offset variable (which acts as a constant on the right hand side of the regression equation.)

In R, with this same dataset from http://data.princeton.edu/wws509/datasets/#smoking (which I'm storing as raw.dat, with variable names age, smoke, pop and dead) this would be:

# Make a factorial version of the smoking variable
raw.dat$smoke.f <- factor(raw.dat$smoke)

# For this variable 1 = non-smoker, 2 = cigar only,
# 3 = cigarettes + cigars, 4 = cigarettes only.

# Regression model:
poisson.smoke <- glm(dead ~ age + smoke.f, offset=log(pop),
family="poisson", data = raw.dat)
summary(poisson.smoke)

# snipped output from R:
# Coefficients:
#              Estimate Std. Error z value Pr(>|z|)
# (Intercept) -3.738877   0.050009 -74.764  < 2e-16 ***
# age          0.333006   0.005591  59.559  < 2e-16 ***
# smoke.f2     0.032927   0.046894   0.702    0.483
# smoke.f3     0.236353   0.038597   6.124 9.15e-10 ***
# smoke.f4     0.437946   0.039803  11.003  < 2e-16 ***


So the coefficient for smoke.f2 (cigar only) compared to non-smokers is now 0.333, exponentiating this gives a rate ratio of 1.033 (95% CI 0.94, 1.13) so no particular evidence in this model of increased risk of lung cancer in cigar-only smokers compared to non-smokers.

Secondly: regarding the lack of interaction term between cigar and cigarette in your original model (and note that including an interaction term is at one level equivalent to treating cigar + cigarette as a distinct group, as per your second model and the example above, which leads to a model with the same number of parameters and the same solution, but some different interpretations of coefficients.)

This lack of interaction might be causing some issues, but with correct specification of the offset variable the coefficient is at least pointing in the correct direction:

# Make two indicator variables for cigar and smoking
raw.dat$cigar <- ifelse(raw.dat$smoking==2 | raw.dat$smoking==3, 1, 0) raw.dat$cigarette <- ifelse(raw.dat$smoking==3 | raw.dat$smoking==4, 1, 0)

# Model treating these as independent:
poisson.sep.smoke <- glm(dead ~ age + cigar + cigarette,  offset=log(pop),
family="poisson", data = raw.dat)

# snipped output from R:
# Coefficients:
#              Estimate Std. Error z value Pr(>|z|)
# (Intercept) -3.648010   0.044886 -81.273  < 2e-16 ***
# age          0.334359   0.005576  59.959  < 2e-16 ***
# cigar       -0.153685   0.021252  -7.232 4.77e-13 ***
# cigarette    0.312363   0.027301  11.442  < 2e-16 ***


And then including an interaction between the two: solution-wise, this is equivalent to having all four levels (with three dummy variables) as per my first model.

poisson.int.smoke <- glm(dead ~ age + cigar * cigarette,  offset=log(pop),
family="poisson", data = raw.dat)
summary(poisson.int.smoke)

# snipped output from R:
# Coefficients:
#                  Estimate Std. Error z value Pr(>|z|)
# (Intercept)     -3.738877   0.050009 -74.764  < 2e-16 ***
# age              0.333006   0.005591  59.559  < 2e-16 ***
# cigar            0.032927   0.046894   0.702    0.483
# cigarette        0.437946   0.039803  11.003  < 2e-16 ***
# cigar:cigarette -0.234519   0.052428  -4.473 7.71e-06 ***


So you can see that the interaction term (cigar:cigarette) is significant: in other words, the effects of cigar smoking and cigarette smoking are dependent on each other. This is the main advantage of specifying an interaction-term based model over the all-levels model, in that the interaction phenomenon (of independence between cigar and cigarette impact) is "automatically" tested in this scenario.

The cigar and cigarette coefficients here correspond to smoke.f2 and smoke.f4 before, and mathematically, the old coefficient for smoke.f3 is the sum of cigar, cigarette and cigar:cigarette terms in this model.

• (+1) Thanks for taking the time to work up all these examples and post them. I think you've covered all the bases. It's curious that group 3 has lower incidence than group 4. Maybe if you smoke a pipe/cigar you tend to be more recreational about your cigarette use than if you smoke only cigarettes. – cardinal Jan 31 '13 at 23:20
• I think you're suspicion is right, it's highly likely that the "mixed" group are regular cigar/pipe smokers (I work with people here who look at tobacco control issues, although I haven't dealt with much data on it myself,) so that would explain the intermediate risk profile for that group. – James Stanley Feb 1 '13 at 0:03
• great! thanks for the assistance James, ill make the corrections to my model. – Cam.Davidson.Pilon Feb 1 '13 at 2:21
• You're welcome -- as I said above, I counted this as good practice for my upcoming teaching :-) – James Stanley Feb 1 '13 at 2:59
• I got the same results. You students are lucky – Cam.Davidson.Pilon Feb 1 '13 at 13:13