Documentation Example 3 for PROC MCMC
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: MCMCEX3 */
/* TITLE: Documentation Example 3 for PROC MCMC */
/* Logistic Regression Model with Diffuse Prior */
/* PRODUCT: STAT */
/* SYSTEM: ALL */
/* KEYS: */
/* PROCS: MCMC */
/* DATA: */
/* */
/* SUPPORT: Fang Chen */
/* REF: PROC MCMC, EXAMPLE 3 */
/* MISC: */
/****************************************************************/
/****************************************************************/
/* Logistic Regression Model with a Diffuse Prior */
/****************************************************************/
title 'Logistic Regression Model with a Diffuse Prior';
data beetles;
input n y x @@;
datalines;
6 0 25.7 8 2 35.9 5 2 32.9 7 7 50.4 6 0 28.3
7 2 32.3 5 1 33.2 8 3 40.9 6 0 36.5 6 1 36.5
6 6 49.6 6 3 39.8 6 4 43.6 6 1 34.1 7 1 37.4
8 2 35.2 6 6 51.3 5 3 42.5 7 0 31.3 3 2 40.6
;
ods graphics on;
proc mcmc data=beetles ntu=1000 nmc=20000 nthin=2 propcov=quanew
diag=(mcse ess) outpost=beetleout seed=246810;
ods select PostSummaries PostIntervals mcse ess TADpanel;
parms (alpha beta) 0;
prior alpha beta ~ normal(0, var = 10000);
p = logistic(alpha + beta*x);
model y ~ binomial(n,p);
run;
proc format;
value betafmt low-0 = 'beta <= 0' 0<-high = 'beta > 0';
run;
proc freq data=beetleout;
tables beta /nocum;
format beta betafmt.;
run;
proc mcmc data=beetles ntu=1000 nmc=20000 nthin=2 propcov=quanew
outpost=beetleout seed=246810 plot=density
monitor=(pi30 ld05 ld50 ld95);
ods select PostSummaries PostIntervals densitypanel;
parms (alpha beta) 0;
begincnst;
c1 = log(0.05 / 0.95);
c2 = -c1;
endcnst;
beginnodata;
prior alpha beta ~ normal(0, var = 10000);
pi30 = logistic(alpha + beta*30);
ld05 = (c1 - alpha) / beta;
ld50 = - alpha / beta;
ld95 = (c2 - alpha) / beta;
endnodata;
pi = logistic(alpha + beta*x);
model y ~ binomial(n,pi);
run;
ods graphics off;
data transout;
set beetleout;
pi30 = logistic(alpha + beta*30);
ld05 = (log(0.05 / 0.95) - alpha) / beta;
ld50 = (log(0.50 / 0.50) - alpha) / beta;
ld95 = (log(0.95 / 0.05) - alpha) / beta;
run;