~john-c/5421/zeroinflatedpoisson.notes
Below are two programs and their output which analyze parameters for
zero-inflated Poisson data. Here the assumption is that the data are
from a mixture of two populations: one which yields only 0's as measurements,
and the other where the measurements are Poisson-distributed counts with
mean lambda.
It is assumed that the mixing parameter is alpha, that is, that a
proportion alpha of the measurements come from Population 1 and a
proportion (1 - alpha) come from Population 2.
Note that people who have measurement 0 can come from either of
the two populations: 0 is a possible value for a Poisson distribution.
Thus the likelihood for an observation Y is either:
alpha + (1 - alpha) * exp(-lambda) if Y = 0, or
(1 - alpha) * lambda^Y * exp(-lambda) / Y! if Y > 0.
The first program is just a Gauss-Newton maximum likelihood program.
The second program is based on the EM algorithm. This program does
not operate on the original data, but on the summary counts. The
EM algorithm is usually applied when some data are missing. A crude
version of the EM idea is to (1) impute the expected values of the missing
data from an initial set of parameters, to yield 'pseudo-complete' data,
and (2) compute the next iteration of the parameters by maximizing the
function resulting from Step (1). These are the 'E' and 'M' steps
respectively.
In this EM example, the missing data are the numbers of people from Population 1
and Population 2 who have measured values of 0. These numbers are denoted
by xa and xb respectively in the program. Note that the data give the total
number of people with the measured value of 0, but they do not tell you
how many came from each population. That is, xa + xb = totalzeroes is
known, but xa and xb themselves are not known. In the E-step, each of
xa and xb is estimated; in the M-step, parameters alpha and lambda are
chosen to maximize the pseudo-complete likelihood.
Both of these program give the same estimates for the parameters
alpha and lambda. They also give the same estimates as can be obtained
from PROC GENMOD using an option on the MODEL statement called: dist = zip.
Note that the Gauss-Newton program converged in 11 iterations, while 27
iterations were required for the EM algorithm. Typically mixture problems
present computational difficulties which often prevent convergence of
Gauss-Newton procedures. The EM algorithm typically converges but often
very slowly.
----------------------------------------------------------------------------
options linesize = 100 ;
footnote "~john-c/5421/zeroinflatedpoisson.sas &sysdate &systime" ;
* ;
* This program estimates parameters for zero-inflated Poisson using ;
* a Gauss-Newton algorithm ... ;
* ;
%let numobs = 300 ;
%let truealpha = .4 ;
%let truelambda = 2.0 ;
%let seed = 20111118 ;
data zeroinfpoisson ;
alpha = &truealpha ;
n = &numobs;
lambda = &truelambda ;
do i = 1 to n ;
r = ranuni(&seed) ;
if r < alpha then do ;
count = 0 ;
output ;
end ;
if r >= alpha then do ;
count = ranpoi(&seed, lambda) ;
output ;
end ;
end ;
run ;
proc freq data = zeroinfpoisson ;
tables count ;
run ;
proc iml ;
title1 'PROC IML MLE for zero-inflated Poisson ...' ;
use zeroinfpoisson ;
read all var {count} into t ;
n = nrow(t) ;
p = 2 ;
eps = 1e-8 ;
diff = 99999 ;
oldl = -1e20 ;
evals = 0 ;
beta = {.5, 2.0} ;
e = j(n, 1, 1) ;
start loglike(n, beta, p, t, l, d1l) ;
l = 0 ;
d1l = j(n, p, 0) ;
do i = 1 to n ;
ti = t[i] ;
alpha = beta[1] ;
lambda = beta[2] ;
emlambda = exp(-lambda) ;
if ti = 0 then like = alpha + (1 - alpha) * exp(-lambda) ;
if ti > 0 then like = (1 - alpha) * (lambda**ti) * exp(-lambda) ;
if ti = 0 then loglike = log(alpha + (1 - alpha) * exp(-lambda)) ;
if ti > 0 then loglike = log(1 - alpha) + ti * log(lambda) - lambda ;
l = l + loglike ;
if ti = 0 then d1ldalpha = (1 - exp(-lambda)) / (alpha + (1 - alpha) * exp(-lambda)) ;
if ti > 0 then d1ldalpha = -1 / (1 - alpha) ;
if ti = 0 then d1ldlambda = -(1 - alpha) * exp(-lambda) / like ;
if ti > 0 then d1ldlambda = ti / lambda - 1 ;
d1l[i, 1] = d1ldalpha ;
d1l[i, 2] = d1ldlambda ;
end ;
finish loglike ;
run loglike(n, beta, p, t, l, d1l) ;
print "Start values:" beta l ;
maxderiv = 99999 ;
do iter = 1 to 50 while (diff > eps & maxderiv > eps) ;
run loglike(n, beta, p, t, l, d1l) ;
delta = inv(d1l` * d1l) * (d1l` * e) ;
oldbeta = beta ;
beta = beta + delta ;
diff = max(abs(delta)) ;
maxderiv = max(abs(d1l` * e)) ;
print iter oldbeta beta l diff maxderiv ;
end ;
if (diff > eps) then do ;
print "No convergence after " iter " iterations ..." ;
end ;
if (diff < eps) then do ;
covar = inv(d1l` * d1l) ;
serrbeta1 = sqrt(covar[1, 1]) ;
serrbeta2 = sqrt(covar[2, 2]) ;
print "---------------------------------------------" ;
print "Number of observations: " &numobs ;
print "Convergence achieved ..." ;
print "Final Parameter Estimates: " ;
beta1hat = beta[1] ; quotient1 = beta1hat / serrbeta1 ;
beta2hat = beta[2] ; quotient2 = beta2hat / serrbeta2 ;
print "beta1hat = " beta1hat " serr(beta1hat) = " serrbeta1 " beta1hat/serrbeta1 = " quotient1 ;
print "beta2hat = " beta2hat " serr(beta2hat) = " serrbeta2 " beta2hat/serrbeta2 = " quotient2 ;
print "True values of probzero and lambda = " &truealpha &truelambda ;
m2logl = -2 * l ;
print "-2 * log(likelihood) = " m2logl ;
print "---------------------------------------------" ;
end ;
quit ;
proc freq data = zeroinfpoisson ;
tables count ;
run ;
proc genmod data = zeroinfpoisson ;
model count = / dist = zip ;
zeromodel / link = logit ;
title 'PROC GENMOD for zero-inflated Poisson ...' ;
run ;
PROC IML MLE for zero-inflated Poisson ... 2
18:44 Monday, November 21, 2011
beta l
Start values: 0.5 -237.5108
2
iter oldbeta beta l diff maxderiv
1 0.5 0.4677629 -237.5108 0.0773609 36.28987
2 2.0773609
iter oldbeta beta l diff maxderiv
2 0.4677629 0.4661239 -236.6476 0.0050742 1.1578301
2.0773609 2.0722867
iter oldbeta beta l diff maxderiv
3 0.4661239 0.4661958 -236.6464 0.000997 0.0455482
2.0722867 2.0732837
iter oldbeta beta l diff maxderiv
4 0.4661958 0.4661806 -236.6463 0.0001956 0.0088474
2.0732837 2.0730881
iter oldbeta beta l diff maxderiv
5 0.4661806 0.4661836 -236.6463 0.0000385 0.0017428
2.0730881 2.0731266
iter oldbeta beta l diff maxderiv
6 0.4661836 0.466183 -236.6463 7.5617E-6 0.0003425
2.0731266 2.073119
iter oldbeta beta l diff maxderiv
7 0.466183 0.4661831 -236.6463 1.4867E-6 0.0000673
2.073119 2.0731205
iter oldbeta beta l diff maxderiv
8 0.4661831 0.4661831 -236.6463 2.9229E-7 0.0000132
2.0731205 2.0731202
~john-c/5421/zeroinflatedpoisson.sas 21NOV11 18:44
PROC IML MLE for zero-inflated Poisson ... 3
18:44 Monday, November 21, 2011
iter oldbeta beta l diff maxderiv
9 0.4661831 0.4661831 -236.6463 5.7465E-8 2.6031E-6
2.0731202 2.0731203
iter oldbeta beta l diff maxderiv
10 0.4661831 0.4661831 -236.6463 1.1298E-8 5.1177E-7
2.0731203 2.0731203
iter oldbeta beta l diff maxderiv
11 0.4661831 0.4661831 -236.6463 2.2212E-9 1.0062E-7
2.0731203 2.0731203
---------------------------------------------
Number of observations: 300
Convergence achieved ...
Final Parameter Estimates:
beta1hat serrbeta1 quotient1
beta1hat = 0.4661831 serr(beta1hat) = 0.0348685 beta1hat/serrbeta1 = 13.369734
beta2hat serrbeta2 quotient2
beta2hat = 2.0731203 serr(beta2hat) = 0.148579 beta2hat/serrbeta2 = 13.952982
True values of probzero and lambda = 0.4 2
m2logl
-2 * log(likelihood) = 473.29268
---------------------------------------------
~john-c/5421/zeroinflatedpoisson.sas 21NOV11 18:44
PROC IML MLE for zero-inflated Poisson ... 4
18:44 Monday, November 21, 2011
The FREQ Procedure
Cumulative Cumulative
count Frequency Percent Frequency Percent
----------------------------------------------------------
0 160 53.33 160 53.33
1 37 12.33 197 65.67
2 46 15.33 243 81.00
3 36 12.00 279 93.00
4 11 3.67 290 96.67
5 9 3.00 299 99.67
6 1 0.33 300 100.00
~john-c/5421/zeroinflatedpoisson.sas 21NOV11 18:44
PROC GENMOD for zero-inflated Poisson ... 5
18:44 Monday, November 21, 2011
The GENMOD Procedure
Model Information
Data Set WORK.ZEROINFPOISSON
Distribution Zero Inflated Poisson
Link Function Log
Dependent Variable count
Number of Observations Read 300
Number of Observations Used 300
Criteria For Assessing Goodness Of Fit
Criterion DF Value Value/DF
Deviance 835.3194
Scaled Deviance 835.3194
Pearson Chi-Square 298 282.4116 0.9477
Scaled Pearson X2 298 282.4116 0.9477
Log Likelihood -236.6463
Full Log Likelihood -417.6597
AIC (smaller is better) 839.3194
AICC (smaller is better) 839.3598
BIC (smaller is better) 846.7270
Algorithm converged.
Analysis Of Maximum Likelihood Parameter Estimates
Standard Wald 95% Confidence Wald
Parameter DF Estimate Error Limits Chi-Square Pr > ChiSq
Intercept 1 0.7291 0.0655 0.6006 0.8575 123.82 <.0001
Scale 0 1.0000 0.0000 1.0000 1.0000
NOTE: The scale parameter was held fixed.
Analysis Of Maximum Likelihood Zero Inflation Parameter Estimates
Standard Wald 95% Confidence Wald
Parameter DF Estimate Error Limits Chi-Square Pr > ChiSq
Intercept 1 -0.1355 0.1389 -0.4077 0.1367 0.95 0.3293
~john-c/5421/zeroinflatedpoisson.sas 21NOV11 18:44
-------------------------------------------------------------------------------------------------
options linesize = 150;
footnote "~john-c/5421/sasprog.sas &sysdate &systime" ;
* This is an EM-algorithm approach to estimating parameters ;
* for zero-inflated Poisson data. ;
proc iml ;
counts = {160, 37, 46, 36, 11, 9, 1} ;
values = {0, 1, 2, 3, 4, 5, 6} ;
p = .5 ;
lambda = 0.8 ;
print "Initial estimates of p, lambda: " p lambda ;
y0 = counts[1] ;
do i = 1 to 50 ;
xa = counts[1] * p / (p + (1 - p) * exp(-lambda)) ;
xb = counts[1] - xa ;
print "Step " i " E-step: XA, XB: " xa xb ;
numvalues = max(values) ;
sumky = 0 ;
sumy = xb ;
loglikeobs = counts[1] * log(p + (1 - p)*exp(-lambda)) ;
loglikecomplete = xa * log(p) + xb * (log(1 - p) - lambda) ;
do j = 2 to numvalues + 1 ;
sumky = sumky + values[j] * counts[j] ;
sumy = sumy + counts[j] ;
loglikeobs = loglikeobs + counts[j] * (log(1 - p) + (j - 1) * log(lambda) - lambda) ;
loglikecomplete = loglikecomplete + counts[j] * (log(1 - p) + (j - 1) * log(lambda) - lambda) ;
end ;
p = xa / (xa + sumy) ;
lambda = sumky / sumy ;
print "Step " i " M-step: p, lambda, loglikeobs loglikecomplete: " p lambda loglikeobs loglikecomplete ;
end ;
quit ;
The SAS System 17:47 Monday, November 21, 2011 1
P LAMBDA
Initial estimates of p, lambda: 0.5 0.8
I XA XB
Step 1 E-step: XA, XB: 110.39592 49.604083
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 1 M-step: p, lambda, loglikeobs loglikecomplete: 0.3679864 1.7510171 -334.6517 -433.7111
I XA XB
Step 2 E-step: XA, XB: 123.25214 36.74786
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 2 M-step: p, lambda, loglikeobs loglikecomplete: 0.4108405 1.878382 -241.5981 -327.8193
I XA XB
Step 3 E-step: XA, XB: 131.23641 28.763591
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 3 M-step: p, lambda, loglikeobs loglikecomplete: 0.4374547 1.967249 -238.3656 -313.7333
I XA XB
Step 4 E-step: XA, XB: 135.61291 24.387095
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 4 M-step: p, lambda, loglikeobs loglikecomplete: 0.452043 2.0196233 -237.1374 -305.4387
I XA XB
Step 5 E-step: XA, XB: 137.82853 22.171474
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 5 M-step: p, lambda, loglikeobs loglikecomplete: 0.4594284 2.0472158 -236.7693 -301.1473
~john-c/5421/sasprog.sas 21NOV11 17:47
The SAS System 17:47 Monday, November 21, 2011 2
I XA XB
Step 6 E-step: XA, XB: 138.90166 21.09834
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 6 M-step: p, lambda, loglikeobs loglikecomplete: 0.4630055 2.060853 -236.6749 -299.0614
I XA XB
Step 7 E-step: XA, XB: 139.40983 20.59017
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 7 M-step: p, lambda, loglikeobs loglikecomplete: 0.4646994 2.0673744 -236.6527 -298.0745
I XA XB
Step 8 E-step: XA, XB: 139.64784 20.35216
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 8 M-step: p, lambda, loglikeobs loglikecomplete: 0.4654928 2.070443 -236.6477 -297.6127
I XA XB
Step 9 E-step: XA, XB: 139.75874 20.241263
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 9 M-step: p, lambda, loglikeobs loglikecomplete: 0.4658625 2.0718758 -236.6466 -297.3977
I XA XB
Step 10 E-step: XA, XB: 139.81028 20.189718
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 10 M-step: p, lambda, loglikeobs loglikecomplete: 0.4660343 2.0725425 -236.6464 -297.2978
I XA XB
Step 11 E-step: XA, XB: 139.83421 20.165788
~john-c/5421/sasprog.sas 21NOV11 17:47
The SAS System 17:47 Monday, November 21, 2011 3
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 11 M-step: p, lambda, loglikeobs loglikecomplete: 0.466114 2.0728522 -236.6464 -297.2514
I XA XB
Step 12 E-step: XA, XB: 139.84532 20.154684
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 12 M-step: p, lambda, loglikeobs loglikecomplete: 0.4661511 2.0729959 -236.6463 -297.2299
I XA XB
Step 13 E-step: XA, XB: 139.85047 20.149532
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 13 M-step: p, lambda, loglikeobs loglikecomplete: 0.4661682 2.0730626 -236.6463 -297.2199
I XA XB
Step 14 E-step: XA, XB: 139.85286 20.147143
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 14 M-step: p, lambda, loglikeobs loglikecomplete: 0.4661762 2.0730935 -236.6463 -297.2153
I XA XB
Step 15 E-step: XA, XB: 139.85397 20.146034
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 15 M-step: p, lambda, loglikeobs loglikecomplete: 0.4661799 2.0731078 -236.6463 -297.2131
I XA XB
Step 16 E-step: XA, XB: 139.85448 20.14552
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 16 M-step: p, lambda, loglikeobs loglikecomplete: 0.4661816 2.0731145 -236.6463 -297.2121
~john-c/5421/sasprog.sas 21NOV11 17:47
The SAS System 17:47 Monday, November 21, 2011 4
I XA XB
Step 17 E-step: XA, XB: 139.85472 20.145282
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 17 M-step: p, lambda, loglikeobs loglikecomplete: 0.4661824 2.0731176 -236.6463 -297.2117
I XA XB
Step 18 E-step: XA, XB: 139.85483 20.145171
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 18 M-step: p, lambda, loglikeobs loglikecomplete: 0.4661828 2.073119 -236.6463 -297.2115
I XA XB
Step 19 E-step: XA, XB: 139.85488 20.14512
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 19 M-step: p, lambda, loglikeobs loglikecomplete: 0.4661829 2.0731197 -236.6463 -297.2114
I XA XB
Step 20 E-step: XA, XB: 139.8549 20.145096
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 20 M-step: p, lambda, loglikeobs loglikecomplete: 0.466183 2.07312 -236.6463 -297.2113
I XA XB
Step 21 E-step: XA, XB: 139.85491 20.145085
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 21 M-step: p, lambda, loglikeobs loglikecomplete: 0.466183 2.0731201 -236.6463 -297.2113
I XA XB
Step 22 E-step: XA, XB: 139.85492 20.14508
~john-c/5421/sasprog.sas 21NOV11 17:47
The SAS System 17:47 Monday, November 21, 2011 5
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 22 M-step: p, lambda, loglikeobs loglikecomplete: 0.4661831 2.0731202 -236.6463 -297.2113
I XA XB
Step 23 E-step: XA, XB: 139.85492 20.145078
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 23 M-step: p, lambda, loglikeobs loglikecomplete: 0.4661831 2.0731202 -236.6463 -297.2113
I XA XB
Step 24 E-step: XA, XB: 139.85492 20.145077
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 24 M-step: p, lambda, loglikeobs loglikecomplete: 0.4661831 2.0731202 -236.6463 -297.2113
I XA XB
Step 25 E-step: XA, XB: 139.85492 20.145076
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 25 M-step: p, lambda, loglikeobs loglikecomplete: 0.4661831 2.0731202 -236.6463 -297.2113
I XA XB
Step 26 E-step: XA, XB: 139.85492 20.145076
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 26 M-step: p, lambda, loglikeobs loglikecomplete: 0.4661831 2.0731203 -236.6463 -297.2113
I XA XB
Step 27 E-step: XA, XB: 139.85492 20.145076
I P LAMBDA LOGLIKEOBS LOGLIKECOMPLETE
Step 27 M-step: p, lambda, loglikeobs loglikecomplete: 0.4661831 2.0731203 -236.6463 -297.2113
~john-c/5421/sasprog.sas 21NOV11 17:47
zeroinflatedpoisson.notes Most recent update: November 21, 2011.