/**************************************************************************************************
 * Program: zip.sas
 * Description: SAS macro to fit a Zero-inflated Poisson (ZIP) Model in nlmixed.
 *              ZIP model has a zero with probability p and Poisson with probability 1-p; 
 *              the Poisson part has mean lambda
 * Macro Parameters:
	ds - input dataset (default=_last_)
	y - dependent/response variable with non-negative integer values
   	xlist - optional indpedent/predictor variable list, separated by space(s)
   	type -  how independent/predictor variables affect the dependent/response variable:
    	    prob = probability of a zero (ie not Poisson)
            mean = mean (lambda) of the Poisson portion
            both = both probability and mean
   	
 * Note: parameter names are generated as follows:
   bp_varname is the parameter name for a variable with varname that affects probility,
   bl_varname is the parameter name for a variable with varname that affects mean;
   these names can be very long depending on how you name your variables.
 * Author: Nancy Cheng
 * Date: 23 May 2005
 * reference code and simulation data: see zip-model code by Dale McLarren posted on SAS-L listserv
****************************************************************************************************/
*options symbolgen mprint;
%macro zip(ds=, y=, xlist=, type=both);
proc nlmixed data=&ds; 
    %let i=1;
	eta_prob = bp_0;
	eta_lambda = bl_0;
	%do %while(%scan(&xlist,&i, %str( )) ne );
		%let xi=%scan(&xlist, &i,  %str( ));
		%if &type eq %str(prob) %then %do;
			eta_prob = eta_prob + bp_&xi * &xi ; 
		%end;
		%else %if &type eq %str(mean) %then %do;
			eta_lambda = eta_lambda + bl_&xi * ξ
		%end;
		%else %do; /* type=both */
			eta_prob = eta_prob + bp_&xi * &xi ; 
			eta_lambda = eta_lambda + bl_&xi * ξ
		%end;
		
		%let i=%eval(&i+1);	
	%end;
    p_0 = exp(eta_prob)/(1 + exp(eta_prob));
	lambda = exp(eta_lambda);

	/* likelihood structure for ZIP. Note that we cannot generate */ 
	/* the log likelihood for the zero response values directly. We */
	/* must generate the likelihood and then take the logarithm due */ 
	/* to the additive nature of the zero likelihood. */ 
	if &y=0 then like = p_0 + (1-p_0)*exp(-lambda);
	if &y=0 then loglike = log(like);
	else loglike = log(1-p_0) + &y*log(lambda) -lambda - lgamma(&y+1);

	/* fit the model */ 
	model &y ~ general(loglike);
	estimate 'p_0' p_0;
	estimate 'lambda' lambda;
	predict exp(loglike)  out=predzip;
run;	
%mend;

*simulated data; 
data test;
	* parameters for zero inflation probability;
	bp_0=-.5; bp_1=0.3;
    * parameters for poisson mean;
	bl_0=.4; bl_1=.2;

	seed=123157;
	do i=1 to 1000;
	    * generate variables affecting probability and mean eta;
		call ranbin(seed, 1, .5, gender);
		call rannor(seed, x);

		* probability of zero inflation, a function of gender only;
		eta_prob=bp_0 + bp_1*gender;
		p_0=exp(eta_prob)/(1+ exp(eta_prob));
		call ranbin(seed, 1, p_0, zero_inflate);

		* poisson mean, a function of x only;
		eta_lambda=bl_0 + bl_1*x;
		lambda=exp(eta_lambda);

		*generate response;
		if zero_inflate then y=0;
		else call ranpoi(seed, lambda, y);
		output;
	end;
	keep y gender x;
run;
%zip(ds=test, y=y, xlist=gender x, type=both);