This SAS Code Fragment page shows examples of simple linear and nonlinear models using propc mcmc. Please note that the page only shows diagnostic plots for the first model. Users should genrate and inspect diagnostic plots for all models. Also, the prior distributions shown here are used as examples and are not necessarily being recommended.
/*************************/
/* linear regression */
/*************************/
options nocenter;
ods graphic on;
proc mcmc data='D:datahsbdemo.sas7bdat' outpost=linout nmc=20000
propcov=quanew seed=7 diag=(mcse ess) thin=5;
ods select PostSummaries PostIntervals mcse ess TADPanel;
parms b0 0 b1 0 b2 0 s2e 1;
prior b:~ normal(mean = 0, var = 1e6);
prior s2e ~ igamma(shape = 3/10, scale = 10/3);
mu = b0 + b1*read+b2*female;
model write ~ normal(mu, var=s2e);
run;
Posterior Summaries
Standard Percentiles
Parameter N Mean Deviation 25% 50% 75%
b0 4000 20.1659 2.6710 18.3856 20.1994 22.0065
b1 4000 0.5669 0.0489 0.5333 0.5668 0.5986
b2 4000 5.4809 1.0063 4.8167 5.4975 6.1475
s2e 4000 51.2323 5.2743 47.7175 50.9084 54.5646
Posterior Intervals
Parameter Alpha Equal-Tail Interval HPD Interval
b0 0.050 14.8050 25.2491 14.7609 25.1714
b1 0.050 0.4758 0.6687 0.4718 0.6624
b2 0.050 3.4522 7.4388 3.4310 7.3898
s2e 0.050 41.7383 62.6536 40.9200 61.5277
Monte Carlo Standard Errors
Standard
Parameter MCSE Deviation MCSE/SD
b0 0.0755 2.6710 0.0283
b1 0.00138 0.0489 0.0283
b2 0.0244 1.0063 0.0242
s2e 0.1514 5.2743 0.0287
Effective Sample Sizes
Correlation
Parameter ESS Time Efficiency
b0 1252.1 3.1946 0.3130
b1 1251.1 3.1971 0.3128
b2 1703.8 2.3477 0.4260
s2e 1213.7 3.2957 0.3034
/**********************************/
/* Binary logistic regression */
/**********************************/
proc mcmc data='D:datahsbdemo.sas7bdat' outpost=logout nmc=20000
propcov=quanew ntu=1000 seed=7 diag=(mcse ess);
ods select PostSummaries PostIntervals mcse ess TADPanel;
parms b0 0 b1 0 b2 0;
prior b: ~ normal(0, prec = 1e-6);
xb = b0 + b1*read+b2*female;
prob = exp(xb)/(1+exp(xb));
ll = log((prob**honors)*((1-prob)**(1-honors)));
model honors ~ general(ll);
run;
Posterior Summaries
Standard Percentiles
Parameter N Mean Deviation 25% 50% 75%
b0 20000 -9.9065 1.4430 -10.8509 -9.8475 -8.9028
b1 20000 0.1487 0.0235 0.1324 0.1483 0.1640
b2 20000 1.1769 0.4318 0.8835 1.1736 1.4643
Posterior Intervals
Parameter Alpha Equal-Tail Interval HPD Interval
b0 0.050 -12.8586 -7.2115 -12.8142 -7.1884
b1 0.050 0.1041 0.1975 0.1023 0.1948
b2 0.050 0.3410 2.0566 0.2742 1.9720
Monte Carlo Standard Errors
Standard
Parameter MCSE Deviation MCSE/SD
b0 0.0312 1.4430 0.0216
b1 0.000513 0.0235 0.0219
b2 0.0102 0.4318 0.0237
Effective Sample Sizes
Correlation
Parameter ESS Time Efficiency
b0 2141.7 9.3383 0.1071
b1 2088.9 9.5745 0.1044
b2 1783.7 11.2127 0.0892
/**************************/
/* Poisson regression */
/**************************/
proc mcmc data='D:datahsbdemo.sas7bdat' outpost=logout nmc=10000
propcov=quanew ntu=1000 seed=7 diag=(mcse ess);
ods select PostSummaries PostIntervals mcse ess TADPanel;
parms b0 0 b1 0 b2 0;
prior b: ~ normal(0, prec = 1e-6);
xb = b0 + b1*read+b2*female;
mu = exp(xb);
ll = awards*log(mu) - mu - lgamma(awards+1);
model awards ~ general(ll);
run;
Posterior Summaries
Standard Percentiles
Parameter N Mean Deviation 25% 50% 75%
b0 10000 -3.1144 0.3385 -3.3419 -3.1108 -2.8729
b1 10000 0.0602 0.00546 0.0564 0.0601 0.0639
b2 10000 0.4825 0.1167 0.4028 0.4813 0.5594
Posterior Intervals
Parameter Alpha Equal-Tail Interval HPD Interval
b0 0.050 -3.7932 -2.4840 -3.7777 -2.4800
b1 0.050 0.0500 0.0710 0.0498 0.0707
b2 0.050 0.2462 0.7096 0.2565 0.7140
Monte Carlo Standard Errors
Standard
Parameter MCSE Deviation MCSE/SD
b0 0.0118 0.3385 0.0347
b1 0.000192 0.00546 0.0352
b2 0.00387 0.1167 0.0332
Effective Sample Sizes
Correlation
Parameter ESS Time Efficiency
b0 830.0 12.0487 0.0830
b1 808.7 12.3660 0.0809
b2 907.9 11.0139 0.0908
/*************************************/
/* negative binomial regression */
/*************************************/
proc mcmc data='D:datahsbdemo.sas7bdat' outpost=nbout nmc=10000
propcov=quanew ntu=1000 seed=7 diag=(mcse ess);
ods select PostSummaries PostIntervals mcse ess TADPanel;
parms b0 0 b1 0 b2 0 alpha .1;
prior b: ~ normal(0, prec = 1e-6);
prior alpha ~ uniform(0, 1);
xb = b0 + b1*read+b2*female;
mu = exp(xb);
m = 1/alpha;
ll = lgamma(awards+m)-lgamma(awards+1)-lgamma(m)
+ awards*log(alpha*mu)-(awards+m)*log(1+alpha*mu);
model awards ~ general(ll);
run;
Posterior Summaries
Standard Percentiles
Parameter N Mean Deviation 25% 50% 75%
b0 10000 -3.4123 0.4275 -3.6806 -3.3943 -3.1108
b1 10000 0.0652 0.00705 0.0603 0.0649 0.0696
b2 10000 0.5242 0.1446 0.4253 0.5230 0.6209
alpha 10000 0.2265 0.0960 0.1587 0.2175 0.2843
Posterior Intervals
Parameter Alpha Equal-Tail Interval HPD Interval
b0 0.050 -4.2933 -2.6202 -4.2736 -2.6182
b1 0.050 0.0520 0.0797 0.0520 0.0795
b2 0.050 0.2563 0.8294 0.2434 0.8062
alpha 0.050 0.0688 0.4304 0.0668 0.4257
Posterior Summaries
Standard Percentiles
Parameter N Mean Deviation 25% 50% 75%
b0 10000 -3.4123 0.4275 -3.6806 -3.3943 -3.1108
b1 10000 0.0652 0.00705 0.0603 0.0649 0.0696
b2 10000 0.5242 0.1446 0.4253 0.5230 0.6209
alpha 10000 0.2265 0.0960 0.1587 0.2175 0.2843
Posterior Intervals
Parameter Alpha Equal-Tail Interval HPD Interval
b0 0.050 -4.2933 -2.6202 -4.2736 -2.6182
b1 0.050 0.0520 0.0797 0.0520 0.0795
b2 0.050 0.2563 0.8294 0.2434 0.8062
alpha 0.050 0.0688 0.4304 0.0668 0.4257