Prev: Soundex singed different values for same names?
Next: Datasetp Retain Question
From: Ryan on 11 Mar 2010 15:30
Dear SAS-L,

I've attempted to run a single hurdle poisson model employing the
experimental mcmc procedure. When constructing the code, I borrowed
heavily from an old post from Dale showing how to run such a model in
the nlmixed procedure. Any feedback on the code I constructed below
would be most appreciated.

Thanks,

Ryan

p.s. The example came from here:

http://support.sas.com/rnd/app/examples/stat/BayesZIP/bayeszip.pdf

--

data catch;
input gender $ age count @@;
if gender = 'F' then do;
female = 1; male = 0;
end;
else do;
female = 0; male = 1;
end;
obs = _N_;
datalines;
F 54 18 M 37 0 F 48 12 M 27 0
M 55 0 M 32 0 F 49 12 F 45 11
M 39 0 F 34 1 F 50 0 M 52 4
M 33 0 M 32 0 F 23 1 F 17 0
F 44 5 M 44 0 F 26 0 F 30 0
F 38 0 F 38 0 F 52 18 M 23 1
F 23 0 M 32 0 F 33 3 M 26 0
F 46 8 M 45 5 M 51 10 F 48 5
F 31 2 F 25 1 M 22 0 M 41 0
M 19 0 M 23 0 M 31 1 M 17 0
F 21 0 F 44 7 M 28 0 M 47 3
M 23 0 F 29 3 F 24 0 M 34 1
F 19 0 F 35 2 M 39 0 M 43 6
;

/* ------------------------------------------ */
/* Single Hurdle Poisson Model */
/* using mcmc */
/* ------------------------------------------ */

proc mcmc data=catch seed =1181 nmc=100000 thin=10 propcov=quanew ;

ods select Parameters PostSummaries PostIntervals tadpanel;

/* Starting Values */
parms betaL0 0 betaL1 0 betaL2 0 betaP0 0 betaP1 0 betaP2 0;

/* Priors */
prior beta: ~ normal(0,var=1000);

/* Linear function for zero-probability model */
eta_prob = betaL0 + betaL1*female*age + betaL2*male*age;
p_0 = 1 / (1 + exp(-eta_prob));

/* Linear function for Poisson probability model */
eta_lambda = betaP0 + betaP1*female*age + betaP2*male*age;
lambda=exp(eta_lambda);
pYyP = exp(-lambda)*(lambda**count) / fact(count);
pY0P = exp(-lambda);

/* Log Likelihood Construction */
if count=0 then
loglike = eta_prob - log (1+exp(eta_prob));
else
loglike = log(1-p_0) - lambda + count*eta_lambda - lgamma(count+1)
- log(1-pY0P);

model count ~ general(loglike);
run;

/* ------------------------------------------ */
/* Single Hurdle Poisson Model */
/* using nlmixed */
/*------------------------------------------ */

proc nlmixed data=catch;

/*Starting Values*/
parms betaL0 0 betaL1 0 betaL2 0 betaP0 0 betaP1 0 betaP2 0;

/* Linear function for zero-probability model */
eta_prob = betaL0 + betaL1*female*age + betaL2*male*age;
p_0 = 1 / (1 + exp(-eta_prob));

/* Linear function for Poisson probability model */
eta_lambda = betaP0 + betaP1*female*age + betaP2*male*age;
lambda=exp(eta_lambda);
pYyP = exp(-lambda)*(lambda**count) / fact(count);
pY0P = exp(-lambda);

/*Log Likelihood Construction*/
if count=0 then
loglike = eta_prob - log (1+exp(eta_prob));
else
loglike = log(1-p_0) - lambda + count*eta_lambda - lgamma(count+1)
- log(1-pY0P);

model count ~ general(loglike);
run;
 | 
Pages: 1
Prev: Soundex singed different values for same names?
Next: Datasetp Retain Question