Non-intuitive answer from a Poisson regression

Question

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.

Answer 1

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.