CQLIM Procedure
Example 12.3 Bayesian Analysis
This example uses the CQLIM procedure to process a small data table.
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 data set bayes_ex into a data table in your session that is associated with the mylib libref. The DATA step assumes that your libref is named mylib, but you can substitute any appropriately defined libref.
data mylib.bayes_ex;
set bayes_ex;
run;
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 = mylib.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 12.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 12.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 | ||||||