/* Adapted from code provide by Dale McLerran of */
/* Fred Hutchinson Cancer Research Center */
/* Demonstration of ZIP without covariates   */
data test;
seed=12345679;
/* Generate a ZIP sample with a mixing proportion */
/* p0 of .5 and a Poisson mean lambda of 2 */
p0=.5;
do i=1 to 100;
/* Generate the Bernoulli mixing random variable */
  call ranbin(seed,1,p0,zero_inflate);
  lambda=2;
/* Generate the response */
  if zero_inflate then y=0;
  else call ranpoi(seed,lambda,y);
  output;
end;
keep y lambda;
run;
/* Fit the model in NLMIXED.  Use the same start */
/* values as in the algorithm that created the file */

proc nlmixed data=test;
parms p0=.5 mu0=2;
if y=0 then do;
prob=p0+(1-p0)*exp(-mu0);
loglike=log(prob); end;
else loglike=log(1-p0)+y*log(mu0)-mu0-lgamma(y+1);
model y~general(loglike);
run;
/* Demonstration of ZIP with a covariate */
/* Fit a mixed  Poisson model in NLMIXED to the Horseshoe */
/* Crab data.  Use start values based on observed proportion of */
/* 0's and slope parameter estimates from the unmixed model. */
/* The commands below assume you have read the horseshoe */
/* crab data into a file called crab in the WORK directory  */
proc freq data=crab;
table satell;
run;
proc nlmixed data=crab;
parms p0=.5 b0=-.5 b1=.5;
mu0=exp(b0+b1*weight);
if satell=0 then do;
prob=p0+(1-p0)*exp(-mu0);
loglike=log(prob);
end;
else loglike=log(1-p0)+satell*log(mu0)-mu0-lgamma(satell+1);
model satell~general(loglike);
run;
/* Compare this model to a default Poisson loglinear model */
proc nlmixed data=crab;
parms b0=-.5 b1=.5;
mu0=exp(b0+b1*weight);
loglike=satell*log(mu0)-mu0-lgamma(satell+1);
model satell~general(loglike);
run;