 |
|
|
|
| |||
|
|
|||||
This section illustrates the use of Poisson models for counts.
We will be using Stata's poisson command, often
followed by estat gof to compute the deviance.
An alternative way to fit these models is to use the
glm command to fit generalized linear models
with Poisson errors and link log, which reports the deviance
and Pearson's chi-squared statistic.
We will use the data from Fiji on children ever born
that appear on Table 4.1 of the lecture notes.
The data are available on our datasets page at
http://data.princeton.edu/wws509/datasets
in a file called ceb.raw.
The file has 70 lines, one for each cell in
the table. Each line has the cell number, numeric codes for
marriage duration, residence and education, the mean and variance
of children ever born and the number of women in the cell.
We will read the data directly from the website and then label the categories of duration, residence and education.
. infile cell dur res educ mean var women using /// > http://data.princeton.edu/wws509/datasets/ceb.raw (70 observations read) . label define dur 1 "0-4" 2 "5-9" 3 "10-14" 4 "15-19" 5 "20-24" 6 "25-29" . label values dur dur . label define res 1 "Suva" 2 "Urban" 3 "Rural" . label values res res . label define educ 1 "none" 2 "lower pri" 3 "upper pri" 4 "sec +" . label values educ educ
Let us do Figure 4.1, plotting the cell variances versus the cell means in a log-log scale, for cells with at least 20 cases: Note that Stata has an option to use log scales, so we don't need to take logs
. graph twoway scatter var mean if women > 20, ///
> xscale(log) yscale(log) ///
> title("Figure 4.1: Mean-Variance Relationship for CEB Data")
. graph export fig41.png, replace
(note: file fig41.png not found)
(file fig41.png written in PNG format)
Clearly the variance increases with the mean. To proceed further we need to create a few additional variables. First we compute the dependent variable, the total number of children born to women in a cell of the table, and the offset, the log of the number of women:
. gen y = round( mean * women, 1) . gen os = log(women)
We will also need dummy variables for the predictors. I also create a few macros to save some typing later on
. // duration . gen dur0509 = dur == 2 . gen dur1014 = dur == 3 . gen dur1519 = dur == 4 . gen dur2024 = dur == 5 . gen dur2529 = dur == 6 . local dur dur0509 dur1014 dur1519 dur2024 dur2529 . // residence . gen urban = res == 2 . gen rural = res == 3 . local res urban rural . // education . gen lowerPri = educ == 2 . gen upperPri = educ == 3 . gen secPlus = educ == 4 . local educ lowerPri upperPri secPlus
We will start by fitting the null model.
. poisson y, offset(os)
Iteration 0: log likelihood = -2080.664
Iteration 1: log likelihood = -2080.664
Poisson regression Number of obs = 70
LR chi2(0) = -0.00
Prob > chi2 = .
Log likelihood = -2080.664 Pseudo R2 = -0.0000
------------------------------------------------------------------------------
y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_cons | 1.376346 .0097119 141.72 0.000 1.357311 1.395381
os | (offset)
------------------------------------------------------------------------------
. estat gof
Goodness-of-fit chi2 = 3731.851
Prob > chi2(69) = 0.0000
Stata reports as "LR chi2" the model chi-squared statistic,
comparing the model at hand with the /null/ model.
To obtain the deviance statistic, or likelihood ratio chi-squared
statistic comparing the model at hand with the /saturated/
model, you have to use the separate post-estimation command
estat gof.
(Or use glm y, family(poisson) offset(os).)
Next we fit the three one-factor models, starting with residence:
. poisson y `res', offset(os)
Iteration 0: log likelihood = -2051.3779
Iteration 1: log likelihood = -2044.3868
Iteration 2: log likelihood = -2044.3778
Iteration 3: log likelihood = -2044.3778
Poisson regression Number of obs = 70
LR chi2(2) = 72.57
Prob > chi2 = 0.0000
Log likelihood = -2044.3778 Pseudo R2 = 0.0174
------------------------------------------------------------------------------
y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
urban | .1442896 .032448 4.45 0.000 .0806926 .2078866
rural | .2280596 .0278321 8.19 0.000 .1735097 .2826095
_cons | 1.204598 .0249922 48.20 0.000 1.155614 1.253581
os | (offset)
------------------------------------------------------------------------------
. estat gof
Goodness-of-fit chi2 = 3659.279
Prob > chi2(67) = 0.0000
Here is education:
. poisson y `educ', offset(os)
Iteration 0: log likelihood = -1588.3352
Iteration 1: log likelihood = -1545.4751
Iteration 2: log likelihood = -1545.2371
Iteration 3: log likelihood = -1545.2371
Poisson regression Number of obs = 70
LR chi2(3) = 1070.85
Prob > chi2 = 0.0000
Log likelihood = -1545.2371 Pseudo R2 = 0.2573
------------------------------------------------------------------------------
y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
lowerPri | -.2117869 .0216769 -9.77 0.000 -.2542729 -.1693008
upperPri | -.6160532 .0288581 -21.35 0.000 -.6726141 -.5594922
secPlus | -1.224676 .0514108 -23.82 0.000 -1.32544 -1.123913
_cons | 1.647278 .0146932 112.11 0.000 1.61848 1.676076
os | (offset)
------------------------------------------------------------------------------
. estat gof
Goodness-of-fit chi2 = 2660.998
Prob > chi2(66) = 0.0000
And here is duration:
. poisson y `dur', offset(os)
Iteration 0: log likelihood = -315.2481
Iteration 1: log likelihood = -297.80021
Iteration 2: log likelihood = -297.77426
Iteration 3: log likelihood = -297.77426
Poisson regression Number of obs = 70
LR chi2(5) = 3565.78
Prob > chi2 = 0.0000
Log likelihood = -297.77426 Pseudo R2 = 0.8569
------------------------------------------------------------------------------
y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
dur0509 | 1.044886 .0523975 19.94 0.000 .9421893 1.147584
dur1014 | 1.444947 .0502397 28.76 0.000 1.346479 1.543416
dur1519 | 1.706756 .0497474 34.31 0.000 1.609253 1.80426
dur2024 | 1.877474 .0496492 37.81 0.000 1.780164 1.974785
dur2529 | 2.078855 .047507 43.76 0.000 1.985743 2.171967
_cons | -.1036046 .0441511 -2.35 0.019 -.1901391 -.01707
os | (offset)
------------------------------------------------------------------------------
. estat gof
Goodness-of-fit chi2 = 166.072
Prob > chi2(64) = 0.0000
These deviances or 'goodness-of-fit chi-2' statistics should agree with the deviances in Table 4.3 of the notes. (You will notice slight differences due to a small change in rounding procedures. The analysis underlying Table 4.3 calculated the number of children ever born in each subgroup as the product of the mean times the number of cases, and retained a few decimals. This Stata run uses the number of children ever born rounded to the nearest integer. If you omit the rounding you will reproduce the results in the lecture notes exactly.)
Here are the two-factor models. First we do the models with main effects only, starting with residence and education (which doesn't really make sense because it omits marriage duration):
. poisson y `res' `educ', offset(os)
Iteration 0: log likelihood = -1604.8781
Iteration 1: log likelihood = -1538.1812
Iteration 2: log likelihood = -1537.9826
Iteration 3: log likelihood = -1537.9826
Poisson regression Number of obs = 70
LR chi2(5) = 1085.36
Prob > chi2 = 0.0000
Log likelihood = -1537.9826 Pseudo R2 = 0.2608
------------------------------------------------------------------------------
y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
urban | .1221596 .03247 3.76 0.000 .0585196 .1857996
rural | .0603765 .0281966 2.14 0.032 .0051122 .1156408
lowerPri | -.2148242 .021831 -9.84 0.000 -.2576121 -.1720363
upperPri | -.6191827 .0291645 -21.23 0.000 -.6763441 -.5620212
secPlus | -1.224269 .0521705 -23.47 0.000 -1.326521 -1.122017
_cons | 1.584686 .028323 55.95 0.000 1.529174 1.640198
os | (offset)
------------------------------------------------------------------------------
. estat gof
Goodness-of-fit chi2 = 2646.489
Prob > chi2(64) = 0.0000
Next we have duration and residence:
. poisson y `dur' `res', offset(os)
Iteration 0: log likelihood = -313.94224
Iteration 1: log likelihood = -275.20402
Iteration 2: log likelihood = -275.07861
Iteration 3: log likelihood = -275.0786
Poisson regression Number of obs = 70
LR chi2(7) = 3611.17
Prob > chi2 = 0.0000
Log likelihood = -275.0786 Pseudo R2 = 0.8678
------------------------------------------------------------------------------
y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
dur0509 | 1.038611 .0524062 19.82 0.000 .9358969 1.141325
dur1014 | 1.435043 .0502608 28.55 0.000 1.336534 1.533552
dur1519 | 1.698132 .0497642 34.12 0.000 1.600596 1.795668
dur2024 | 1.867528 .0496742 37.60 0.000 1.770168 1.964887
dur2529 | 2.073802 .047513 43.65 0.000 1.980678 2.166925
urban | .1189459 .032464 3.66 0.000 .0553176 .1825742
rural | .18192 .0278516 6.53 0.000 .127332 .2365081
_cons | -.2347572 .0493459 -4.76 0.000 -.3314734 -.1380409
os | (offset)
------------------------------------------------------------------------------
. estat gof
Goodness-of-fit chi2 = 120.6806
Prob > chi2(62) = 0.0000
And finally duration and education:
. poisson y `dur' `educ', offset(os)
Iteration 0: log likelihood = -487.16426
Iteration 1: log likelihood = -265.45863
Iteration 2: log likelihood = -264.8343
Iteration 3: log likelihood = -264.83424
Iteration 4: log likelihood = -264.83424
Poisson regression Number of obs = 70
LR chi2(8) = 3631.66
Prob > chi2 = 0.0000
Log likelihood = -264.83424 Pseudo R2 = 0.8727
------------------------------------------------------------------------------
y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
dur0509 | .9964546 .0527378 18.89 0.000 .8930905 1.099819
dur1014 | 1.368604 .0510511 26.81 0.000 1.268545 1.468662
dur1519 | 1.607787 .0511403 31.44 0.000 1.507554 1.70802
dur2024 | 1.778477 .0511775 34.75 0.000 1.678171 1.878783
dur2529 | 1.96267 .0498885 39.34 0.000 1.86489 2.060449
lowerPri | .0105696 .0223741 0.47 0.637 -.0332829 .0544221
upperPri | -.1228681 .0304484 -4.04 0.000 -.1825459 -.0631903
secPlus | -.3602592 .0540842 -6.66 0.000 -.4662624 -.254256
_cons | .0181279 .0485744 0.37 0.709 -.0770762 .1133321
os | (offset)
------------------------------------------------------------------------------
. estat gof
Goodness-of-fit chi2 = 100.1919
Prob > chi2(61) = 0.0012
Now we add two-factor interactions.
To save a few trees I will use quietly to suppress
the detailed output of the poisson command,
and will print only the deviance using estat gof.
Because we are not looking at the parameter estimates
I will use the xi option to have Stata generate
dummy variables for the main effects and interactions. When
you specify an interaction make sure you leave no spaces around
the *. (If I wanted to print the coefficients I would probably
generate my own dummy variables to improve readability.)
. quietly xi: poisson y i.res*i.educ, offset(os)
. estat gof
Goodness-of-fit chi2 = 2625.888
Prob > chi2(58) = 0.0000
. quietly xi: poisson y i.dur*i.res, offset(os)
. estat gof
Goodness-of-fit chi2 = 108.8968
Prob > chi2(52) = 0.0000
. quietly xi: poisson y i.dur*i.educ, offset(os)
. estat gof
Goodness-of-fit chi2 = 84.53072
Prob > chi2(46) = 0.0005
The best fit so far includes duration and education but has significant lack of fit with a chi-squared of 84.5 on 46 d.f.
Ready for the three-factor models? We start with the three-factor additive model. Compare the deviance with the value listed in Table 4.3 and the parameter estimates with the results in Table 4.4.
. poisson y `dur' `res' `educ', offset(os)
Iteration 0: log likelihood = -623.59688
Iteration 1: log likelihood = -252.64903
Iteration 2: log likelihood = -250.07248
Iteration 3: log likelihood = -250.07108
Iteration 4: log likelihood = -250.07108
Poisson regression Number of obs = 70
LR chi2(10) = 3661.19
Prob > chi2 = 0.0000
Log likelihood = -250.07108 Pseudo R2 = 0.8798
------------------------------------------------------------------------------
y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
dur0509 | .9969348 .0527437 18.90 0.000 .8935591 1.100311
dur1014 | 1.369395 .0510688 26.81 0.000 1.269302 1.469488
dur1519 | 1.613757 .0511949 31.52 0.000 1.513417 1.714097
dur2024 | 1.784911 .0512138 34.85 0.000 1.684534 1.885288
dur2529 | 1.976405 .0500341 39.50 0.000 1.87834 2.07447
urban | .1124186 .0324963 3.46 0.001 .048727 .1761102
rural | .1516602 .0283292 5.35 0.000 .096136 .2071845
lowerPri | .0229728 .0226563 1.01 0.311 -.0214327 .0673783
upperPri | -.1012738 .0309871 -3.27 0.001 -.1620073 -.0405402
secPlus | -.3101495 .0552107 -5.62 0.000 -.4183605 -.2019386
_cons | -.1170972 .0549118 -2.13 0.033 -.2247222 -.0094721
os | (offset)
------------------------------------------------------------------------------
. estat gof
Goodness-of-fit chi2 = 70.66559
Prob > chi2(59) = 0.1421
This model appears to fit the data, with a deviance of on d.f. Please refer to the lecture notes for the interpretation of the parameter estimtes.
Let us now check all possible models with interactions.
Make sure you know how to use the following output to test,
for example, whether we need to add a duration by education
interaction. I will use quietly again to save space.
If you are trying these commands interactively you will
probably not want to supress the output.
. quietly xi: poisson y i.dur i.educ*i.res, offset(os)
. estat gof
Goodness-of-fit chi2 = 59.92104
Prob > chi2(53) = 0.2391
. quietly xi: poisson y i.dur*i.res i.educ, offset(os)
. estat gof
Goodness-of-fit chi2 = 57.13525
Prob > chi2(49) = 0.1986
. quietly xi: poisson y i.dur*i.educ i.res, offset(os)
. estat gof
Goodness-of-fit chi2 = 54.80171
Prob > chi2(44) = 0.1274
. quietly xi: poisson y i.dur*i.educ i.res*i.educ, offset(os)
. estat gof
Goodness-of-fit chi2 = 44.52355
Prob > chi2(38) = 0.2163
. quietly xi: poisson y i.dur*i.res i.res*i.educ, offset(os)
. estat gof
Goodness-of-fit chi2 = 44.31134
Prob > chi2(43) = 0.4161
. quietly xi: poisson y i.dur*i.res i.dur*i.educ, offset(os)
. estat gof
Goodness-of-fit chi2 = 42.65186
Prob > chi2(34) = 0.1467
. quietly xi: poisson y i.dur*i.educ i.dur*i.res i.res*i.educ, offset(os)
. estat gof
Goodness-of-fit chi2 = 30.85619
Prob > chi2(28) = 0.3235
It should be clear from the list of deviances that the additive model does a very good job indeed.
Some of these models may fail in older versions of Stata, which
by default allowed up to 40 parameters per model. The solution
is to increase the maximum via set matsize 60. Stata 9
sets the default to 200, which is more than we need for all of
these models.
If you look at the detailed output you may note that Stata claims
to drop some variables due to multicollinearity. This occurs because
the xi prefix is not terribly smart in its handling of factors
involved in more than one interaction, and apparently tries to include
the main effects twice. The variables dropped are the copies and the
originals are indeed included in the model, as you can verify by
inspecting the listing.
Continue with the addendum on
Extra-Poisson Variation (or Over-Dispersion).
Copyright © Germán Rodríguez, 1993-2006.
Please send feedback to grodri@princeton.edu