CQLIM Procedure
Example 12.5 Bayesian Analysis with Multiple Chains
This example shows how to perform Bayesian analysis with multiple chains in the CQLIM procedure by using a small data table.
The following DATA step creates the Christenson Associates airline data, a frequently cited data set (Greene 2000). The data measure the costs, prices of inputs, and utilization rates for six airlines from 1970 to 1984. This example analyzes the log transformations of cost (C), quantity (Q), and price (PF) and the untransformed load factor (LF). These statements assume that your libref is named mylib, but you can substitute any appropriately defined libref.
data mylib.airline;
input Obs AirlineID T C Q PF LF;
Year = T + 1969;
lC = log(C);
lQ = log(Q);
lPF = log(PF);
label lC = "Log Transformation of Costs";
label lQ = "Log Transformation of Quantity";
label lPF = "Log Transformation of Price of Fuel";
label LF = "Load Factor (utilization index)";
datalines;
1 1 1 1140640 0.95276 106650 0.53449
2 1 2 1215690 0.98676 110307 0.53233
3 1 3 1309570 1.09198 110574 0.54774
4 1 4 1511530 1.17578 121974 0.54085
5 1 5 1676730 1.16017 196606 0.59117
... more lines ...
The following statements estimate a simple linear regression model in the Bayesian framework by using the No-U-Turn Sampler (NUTS). NUTS is specified in the SAMPLER= option, and several values controlling the sampler can be set as arguments to SAMPLER=NUTS. The NCHAIN= option controls the number of Markov chains, and the NSAMPLE= option controls the size of the posterior sample per chain. PROC CQLIM automatically aggregates the results of each of the chains to create the posterior statistics and MCMC diagnostics tables that you request. The BYCHAIN suboption of the STATISTICS= and DIAGNOSTICS= options additionally produces those tables separately for each chain. This example produces the "Posterior Summaries" table, the "MCMC Diagnostic Summaries" table, and the "Heidelberger-Welch Diagnostics" table for each chain as well as for all chains combined.
proc cqlim data = mylib.airline;
model lC = lQ lPF LF;
bayes seed = 2358 nsample = 2000 nchain = 4
sampler = nuts(ntune = 2000)
statistics(bychain) = summary
diagnostics(bychain) = (summary heidel);
run;
Output 12.5.1 shows the Bayesian estimation results. The "Posterior Summaries" table contains the estimates of the posterior mean and posterior quantiles computed using all four chains. The other posterior summary statistics tables that are produced by the STATISTICS= option would similarly be computed using all four chains. The "MCMC Diagnostic Summaries" table contains basic convergence diagnostics aggregated across all four chains. The "Autocorrelation Diagnostics" table would also be aggregated for multiple chains. For more information about how this aggregation is performed, see the section Using Multiple Chains in Chapter 2, Introduction to Bayesian Analysis Procedures. The "Heidelberger-Welch Diagnostics" table displays a summary of all the Heidelberger-Welch tests performed for each chain. Specifically, it shows the number of chains for which each test passed for each parameter, as well as the largest number of iterations that were discarded for the stationarity test. The "Geweke Diagnostics" table and the "Raftery-Lewis Diagnostics" table both would be summarized similarly.
Output 12.5.1: Posterior Summaries and MCMC Diagnostics
| Posterior Summaries | ||||||||
|---|---|---|---|---|---|---|---|---|
| Parameter | N | Mean | Standard Deviation | Percentiles | ||||
| 2.5% | 25% | 50% | 75% | 97.5% | ||||
| Intercept | 8000 | 9.5189 | 0.2369 | 9.0576 | 9.3602 | 9.5200 | 9.6780 | 9.9897 |
| lQ | 8000 | 0.8834 | 0.0136 | 0.8565 | 0.8743 | 0.8835 | 0.8924 | 0.9110 |
| lPF | 8000 | 0.4544 | 0.0206 | 0.4146 | 0.4397 | 0.4543 | 0.4685 | 0.4936 |
| LF | 8000 | -1.6387 | 0.3551 | -2.3521 | -1.8720 | -1.6356 | -1.3933 | -0.9782 |
| _Sigma | 8000 | 0.1262 | 0.00994 | 0.1085 | 0.1192 | 0.1255 | 0.1325 | 0.1466 |
| MCMC Diagnostic Summaries | |||||
|---|---|---|---|---|---|
| Parameter | MCSE | MCSE/SD | ESS | Autocorrelation Time | ESS/N |
| Intercept | 0.00634 | 0.0268 | 1394.9 | 5.7350 | 0.1744 |
| lQ | 0.000358 | 0.0263 | 1449.6 | 5.5187 | 0.1812 |
| lPF | 0.000536 | 0.0261 | 1469.3 | 5.4449 | 0.1837 |
| LF | 0.00865 | 0.0244 | 1686.1 | 4.7448 | 0.2108 |
| _Sigma | 0.000216 | 0.0218 | 2111.8 | 3.7882 | 0.2640 |
| Heidelberger-Welch Diagnostics | |||||
|---|---|---|---|---|---|
| Parameter | Stationarity Test | Half-Width Test | |||
| Chains Tested | Tests Passed | Largest Iterations Discarded | Chains Tested | Tests Passed | |
| Intercept | 4 | 4 | 200 | 4 | 4 |
| lQ | 4 | 4 | 0 | 4 | 4 |
| lPF | 4 | 4 | 200 | 4 | 4 |
| LF | 4 | 4 | 0 | 4 | 4 |
| _Sigma | 4 | 4 | 0 | 4 | 4 |
Detailed results for each chain are available in Output 12.5.2 as a result of the BYCHAIN options. The "Posterior Summaries by Chain" table contains the estimates of the posterior mean and posterior quantiles computed separately for each chain. Similarly, the "MCMC Diagnostic Summaries by Chain" table contains basic convergence diagnostics computed for each chain, and the "Heidelberger-Welch Diagnostics by Chain" table contains the Heidelberger-Welch diagnostics computed for each chain. For more information about the differences between the tables in Output 12.5.1 and Output 12.5.2, see the section Using Multiple Chains in Chapter 2, Introduction to Bayesian Analysis Procedures.
Output 12.5.2: Posterior Summaries and MCMC Diagnostics by Chain
| Posterior Summaries by Chain | |||||||||
|---|---|---|---|---|---|---|---|---|---|
| Chain | Parameter | N | Mean | Standard Deviation | Percentiles | ||||
| 2.5% | 25% | 50% | 75% | 97.5% | |||||
| 1 | Intercept | 2000 | 9.5273 | 0.2426 | 9.0619 | 9.3645 | 9.5192 | 9.6952 | 10.0305 |
| lQ | 2000 | 0.8832 | 0.0132 | 0.8576 | 0.8753 | 0.8830 | 0.8919 | 0.9118 | |
| lPF | 2000 | 0.4524 | 0.0209 | 0.4128 | 0.4383 | 0.4525 | 0.4667 | 0.4925 | |
| LF | 2000 | -1.6099 | 0.3489 | -2.3108 | -1.8336 | -1.6070 | -1.3769 | -0.9698 | |
| _Sigma | 2000 | 0.1263 | 0.0107 | 0.1086 | 0.1185 | 0.1254 | 0.1334 | 0.1471 | |
| 2 | Intercept | 2000 | 9.5156 | 0.2507 | 9.0271 | 9.3435 | 9.5209 | 9.6835 | 10.0036 |
| lQ | 2000 | 0.8837 | 0.0147 | 0.8538 | 0.8735 | 0.8839 | 0.8932 | 0.9138 | |
| lPF | 2000 | 0.4545 | 0.0204 | 0.4128 | 0.4401 | 0.4543 | 0.4688 | 0.4917 | |
| LF | 2000 | -1.6347 | 0.3751 | -2.4184 | -1.8903 | -1.6307 | -1.3722 | -0.9120 | |
| _Sigma | 2000 | 0.1266 | 0.00978 | 0.1091 | 0.1195 | 0.1257 | 0.1324 | 0.1476 | |
| 3 | Intercept | 2000 | 9.5163 | 0.2183 | 9.0936 | 9.3677 | 9.5193 | 9.6629 | 9.9285 |
| lQ | 2000 | 0.8829 | 0.0136 | 0.8565 | 0.8731 | 0.8830 | 0.8927 | 0.9092 | |
| lPF | 2000 | 0.4552 | 0.0204 | 0.4168 | 0.4404 | 0.4553 | 0.4691 | 0.4944 | |
| LF | 2000 | -1.6543 | 0.3499 | -2.3549 | -1.8776 | -1.6564 | -1.4140 | -0.9975 | |
| _Sigma | 2000 | 0.1259 | 0.00976 | 0.1077 | 0.1191 | 0.1256 | 0.1325 | 0.1456 | |
| 4 | Intercept | 2000 | 9.5165 | 0.2348 | 9.0438 | 9.3565 | 9.5211 | 9.6735 | 9.9813 |
| lQ | 2000 | 0.8837 | 0.0129 | 0.8594 | 0.8753 | 0.8839 | 0.8917 | 0.9105 | |
| lPF | 2000 | 0.4553 | 0.0204 | 0.4163 | 0.4403 | 0.4550 | 0.4688 | 0.4949 | |
| LF | 2000 | -1.6559 | 0.3440 | -2.3441 | -1.8875 | -1.6461 | -1.4118 | -1.0155 | |
| _Sigma | 2000 | 0.1259 | 0.00951 | 0.1085 | 0.1195 | 0.1254 | 0.1320 | 0.1464 | |
| MCMC Diagnostic Summaries by Chain | ||||||
|---|---|---|---|---|---|---|
| Chain | Parameter | MCSE | MCSE/SD | ESS | Autocorrelation Time | ESS/N |
| 1 | Intercept | 0.0147 | 0.0604 | 274.1 | 7.2954 | 0.1371 |
| lQ | 0.000784 | 0.0591 | 285.9 | 6.9958 | 0.1429 | |
| lPF | 0.00125 | 0.0598 | 279.8 | 7.1470 | 0.1399 | |
| LF | 0.0175 | 0.0502 | 397.1 | 5.0365 | 0.1986 | |
| _Sigma | 0.000527 | 0.0495 | 408.5 | 4.8960 | 0.2042 | |
| 2 | Intercept | 0.0132 | 0.0528 | 359.1 | 5.5694 | 0.1796 |
| lQ | 0.000730 | 0.0498 | 403.5 | 4.9562 | 0.2018 | |
| lPF | 0.00103 | 0.0505 | 392.8 | 5.0921 | 0.1964 | |
| LF | 0.0186 | 0.0495 | 408.5 | 4.8964 | 0.2042 | |
| _Sigma | 0.000434 | 0.0444 | 506.5 | 3.9486 | 0.2533 | |
| 3 | Intercept | 0.0117 | 0.0535 | 349.4 | 5.7249 | 0.1747 |
| lQ | 0.000697 | 0.0512 | 381.4 | 5.2442 | 0.1907 | |
| lPF | 0.000990 | 0.0486 | 423.4 | 4.7240 | 0.2117 | |
| LF | 0.0174 | 0.0498 | 403.6 | 4.9559 | 0.2018 | |
| _Sigma | 0.000380 | 0.0389 | 659.9 | 3.0308 | 0.3300 | |
| 4 | Intercept | 0.0110 | 0.0466 | 459.7 | 4.3503 | 0.2299 |
| lQ | 0.000638 | 0.0494 | 409.9 | 4.8787 | 0.2050 | |
| lPF | 0.00100 | 0.0491 | 415.2 | 4.8167 | 0.2076 | |
| LF | 0.0156 | 0.0452 | 488.9 | 4.0905 | 0.2445 | |
| _Sigma | 0.000385 | 0.0405 | 610.2 | 3.2776 | 0.3051 | |
| Heidelberger-Welch Diagnostics by Chain | |||||||||
|---|---|---|---|---|---|---|---|---|---|
| Chain | Parameter | Stationarity Test | Half-Width Test | ||||||
| Cramer-von Mises Stat | p-Value | Test Outcome | Iterations Discarded | Half-Width | Mean | Relative Half-Width | Test Outcome | ||
| 1 | Intercept | 0.2967 | 0.1381 | Passed | 200 | 0.0315 | 9.5184 | 0.00331 | Passed |
| lQ | 0.0679 | 0.7645 | Passed | 0 | 0.00154 | 0.8832 | 0.00174 | Passed | |
| lPF | 0.3602 | 0.0923 | Passed | 200 | 0.00252 | 0.4533 | 0.00555 | Passed | |
| LF | 0.0853 | 0.6617 | Passed | 0 | 0.0362 | -1.6099 | -0.0225 | Passed | |
| _Sigma | 0.1054 | 0.5592 | Passed | 0 | 0.00103 | 0.1263 | 0.00818 | Passed | |
| 2 | Intercept | 0.0629 | 0.7955 | Passed | 0 | 0.0284 | 9.5156 | 0.00299 | Passed |
| lQ | 0.1841 | 0.3005 | Passed | 0 | 0.00149 | 0.8837 | 0.00168 | Passed | |
| lPF | 0.0904 | 0.6339 | Passed | 0 | 0.00215 | 0.4545 | 0.00473 | Passed | |
| LF | 0.0904 | 0.6339 | Passed | 0 | 0.0395 | -1.6347 | -0.0241 | Passed | |
| _Sigma | 0.1082 | 0.5459 | Passed | 0 | 0.000934 | 0.1266 | 0.00738 | Passed | |
| 3 | Intercept | 0.0896 | 0.6384 | Passed | 0 | 0.0226 | 9.5163 | 0.00238 | Passed |
| lQ | 0.1032 | 0.5692 | Passed | 0 | 0.00143 | 0.8829 | 0.00162 | Passed | |
| lPF | 0.1125 | 0.5272 | Passed | 0 | 0.00193 | 0.4552 | 0.00425 | Passed | |
| LF | 0.0255 | 0.9886 | Passed | 0 | 0.0357 | -1.6543 | -0.0216 | Passed | |
| _Sigma | 0.0666 | 0.7726 | Passed | 0 | 0.000817 | 0.1259 | 0.00649 | Passed | |
| 4 | Intercept | 0.0374 | 0.9461 | Passed | 0 | 0.0228 | 9.5165 | 0.00239 | Passed |
| lQ | 0.0968 | 0.6006 | Passed | 0 | 0.00132 | 0.8837 | 0.00150 | Passed | |
| lPF | 0.0896 | 0.6382 | Passed | 0 | 0.00210 | 0.4553 | 0.00461 | Passed | |
| LF | 0.2274 | 0.2201 | Passed | 0 | 0.0310 | -1.6559 | -0.0187 | Passed | |
| _Sigma | 0.0565 | 0.8357 | Passed | 0 | 0.000657 | 0.1259 | 0.00522 | Passed | |