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
;
ods graphics on;
proc mcmc data=catch seed=1181 nmc=100000 thin=10
propcov=quanew monitor =(_parms_ Pearson);
ods select Parameters PostSummaries PostIntervals tadpanel;
parms beta0 0 beta1 0 beta2 0;
prior beta: ~ normal(0,var=1000);
mu = exp(beta0 + beta1*female*age + beta2*male*age);
model count ~ poisson(mu);
if obs = 1 then Pearson = 0;
Pearson = Pearson + ((count - mu)**2/mu);
run;
ods graphics off;
ods graphics on;
proc mcmc data=catch seed =1181 nmc=100000 thin=10
propcov=quanew monitor =(_parms_ Pearson);
ods select Parameters PostSummaries PostIntervals tadpanel;
parms beta0 0 beta1 0 beta2 0 eta .3;
prior beta: ~ normal(0,var=1000);
prior eta ~ uniform(0,1);
mu=exp(beta0 + beta1*female*age + beta2*male*age);
llike=log(eta*(count eq 0) + (1-eta)*pdf("poisson",count,mu));
model general(llike);
if obs = 1 then Pearson = 0;
mean = (1 - eta)*mu;
sigma2 = (1 - eta)*mu*(1 + eta*mu);
Pearson = Pearson + ((count - mean)**2/sigma2);
run;
ods graphics off;