The CQLIM Procedure
Example 11.3 Bayesian Analysis
This example uses the CQLIM procedure to process a small data table in the distributed computing environment.
The following DATA step generates a data set that contains 1,000 observations from a censored model. The model contains eight variables.
data bayes_ex;
call streaminit(12345);
array vars x1-x7;
array parms{7} (3 4 2 4 -3 -5 -3);
intercept1 = 0;
intercept2 = 4;
intercept3 = 10;
do i = 1 to 1000;
sum_xb = 0;
do j = 1 to 7;
vars[j] = rand('NORMAL', 0, 1);
sum_xb = sum_xb + parms[j] * vars[j];
end;
g = rand('NORMAL', 0, 1);
if g < -0.5 then intercept = intercept1;
if g >= -0.5 then intercept = intercept2;
if g > 0.5 then intercept = intercept3;
y = intercept + sum_xb + 10 * rand('NORMAL', 0, 1);
if y > 400 then y = 400;
if y < 0 then y = 0;
if g < -0.5 then group = 0;
if g >= -0.5 then group = 1;
if g > 0.5 then group = 2;
output;
end;
keep y x1-x7 group;
run;
You can load the bayes_ex data set into your CAS session by naming your CAS engine libref in the first statement of the following DATA step:
data mycas.bayes_ex;
set bayes_ex;
run;
These statements assume that your CAS engine libref is named mycas, but you can substitute any appropriately defined libref.
The following statements estimate a censored model by using Bayesian methods. Note that the BAYES statement controls the settings for the MCMC sampler and determines which Bayesian output is produced. The NCHAIN= option controls the number of Markov chains, and the NSAMPLE= option controls the size of the posterior sample per chain. The PRIOR statements set the prior distributions for the parameters.
proc cqlim data = mycas.bayes_ex;
class group;
model y = x1-x7 group / censored(lb = 0 ub = 400);
bayes priorsummary(shownames) seed = 72342 nsample = 10000 nchain = 4
sampler = rwm(ntune = 100 ntunestage = 20);
prior intercept ~ normal(mean = 0, sd = 100);
prior x1-x5 ~ normal(mean = 0, sd = 10);
prior x6 ~ t(loc = 0, scale = 10, df = 3, lower = 0);
prior group_0 ~ normal(mean = -4, var = 100);
prior _sigma ~ normal(mean = 0, sd = 100);
run;
Output 11.3.1 shows the Bayesian estimation results for the censored model. The "Prior Summary" table shows detailed information about prior distributions for the parameters in the model. The "Posterior Summaries" table contains two point estimates of the parameters from the posterior samples: the posterior mean and the posterior median (50th percentile). The "MCMC Diagnostics Summaries" table contains basic convergence diagnostics to check whether the Markov chain has converged and whether the sample size is sufficient.
Output 11.3.1: Bayesian Estimation of Censored Model
| Prior Summary | ||||||||||
|---|---|---|---|---|---|---|---|---|---|---|
| Parameter | Name | Prior | Bounds | Hyperparameters | ||||||
| Lower | Upper | Name | Value | Name | Value | Name | Value | |||
| Intercept | Intercept | Normal | -Infty | Infty | Mean | 0 | Variance | 10000 | ||
| x1 | x1 | Normal | -Infty | Infty | Mean | 0 | Variance | 100 | ||
| x2 | x2 | Normal | -Infty | Infty | Mean | 0 | Variance | 100 | ||
| x3 | x3 | Normal | -Infty | Infty | Mean | 0 | Variance | 100 | ||
| x4 | x4 | Normal | -Infty | Infty | Mean | 0 | Variance | 100 | ||
| x5 | x5 | Normal | -Infty | Infty | Mean | 0 | Variance | 100 | ||
| x6 | x6 | Truncated t | 0 | Infty | Location | 0 | Scale | 10 | DF | 3 |
| x7 | x7 | Normal | -Infty | Infty | Mean | 0 | Variance | 1000000 | ||
| group 0 | group_0 | Normal | -Infty | Infty | Mean | -4 | Variance | 100 | ||
| group 1 | group_1 | Normal | -Infty | Infty | Mean | 0 | Variance | 1000000 | ||
| _Sigma | _Sigma | Truncated Normal | 0 | Infty | Mean | 0 | Variance | 10000 | ||
| Posterior Summaries | ||||||||
|---|---|---|---|---|---|---|---|---|
| Parameter | N | Mean | Standard Deviation | Percentiles | ||||
| 2.5% | 25% | 50% | 75% | 97.5% | ||||
| Intercept | 40000 | 10.2101 | 0.6881 | 8.8046 | 9.7512 | 10.2263 | 10.6837 | 11.5101 |
| x1 | 40000 | 2.9478 | 0.3737 | 2.2276 | 2.7001 | 2.9466 | 3.1891 | 3.7152 |
| x2 | 40000 | 3.6635 | 0.3549 | 2.9727 | 3.4211 | 3.6657 | 3.9017 | 4.3697 |
| x3 | 40000 | 2.3640 | 0.3546 | 1.7168 | 2.1229 | 2.3415 | 2.5945 | 3.1004 |
| x4 | 40000 | 3.7139 | 0.3440 | 3.0879 | 3.4800 | 3.6988 | 3.9391 | 4.4343 |
| x5 | 40000 | -3.1058 | 0.3693 | -3.8217 | -3.3589 | -3.1012 | -2.8553 | -2.3860 |
| x6 | 40000 | 0.0282 | 0.0281 | 0.000812 | 0.00822 | 0.0197 | 0.0390 | 0.1032 |
| x7 | 40000 | -2.4930 | 0.3636 | -3.1820 | -2.7399 | -2.4991 | -2.2469 | -1.7793 |
| group 0 | 40000 | -9.4894 | 0.9844 | -11.3458 | -10.1500 | -9.5118 | -8.8705 | -7.4648 |
| group 1 | 40000 | -7.4007 | 0.8598 | -9.0624 | -7.9876 | -7.3987 | -6.8169 | -5.6735 |
| _Sigma | 40000 | 10.4933 | 0.2834 | 9.9736 | 10.2973 | 10.4818 | 10.6785 | 11.0791 |
| MCMC Diagnostic Summaries | ||||||
|---|---|---|---|---|---|---|
| Parameter | MCSE | MCSE/SD | ESS | Autocorrelation Time | ESS/N | |
| Intercept | * | 0.0539 | 0.0783 | 163.1 | 245.3 | 0.00408 |
| x1 | * | 0.0240 | 0.0643 | 242.1 | 165.2 | 0.00605 |
| x2 | * | 0.0206 | 0.0580 | 297.1 | 134.6 | 0.00743 |
| x3 | * | 0.0190 | 0.0536 | 348.6 | 114.7 | 0.00872 |
| x4 | * | 0.0215 | 0.0625 | 255.7 | 156.4 | 0.00639 |
| x5 | * | 0.0211 | 0.0571 | 306.3 | 130.6 | 0.00766 |
| x6 | * | 0.00128 | 0.0454 | 485.4 | 82.4062 | 0.0121 |
| x7 | * | 0.0236 | 0.0650 | 236.3 | 169.2 | 0.00591 |
| group 0 | * | 0.0660 | 0.0670 | 222.4 | 179.8 | 0.00556 |
| group 1 | * | 0.0606 | 0.0705 | 201.4 | 198.6 | 0.00503 |
| _Sigma | * | 0.0162 | 0.0571 | 306.4 | 130.6 | 0.00766 |
| *Autocorrelation Remains | ||||||