The CNTSELECT Procedure

This example shows how to use the CNTSELECT procedure to estimate a model in the Bayesian framework by using the No-U-Turn Sampler (NUTS) with multiple Markov chains.

The following DATA step creates the data set DocVisit, which contains information for approximately 5,000 Australian individuals about the number and possible determinants of doctor visits that were made during a two-week interval. This data set contains a subset of variables that are taken from the data set Racd3 that was used by Cameron and Trivedi (1998).

data docvisit;
input sex age agesq income levyplus freepoor
      freerepa illness actdays hscore chcond1 chcond2
      doctorco nondocco hospadmi hospad medicine prescrib nonpresc;
datalines;
1 0.19 0.0361 0.55 1 0 0 1 4 1 0 0 1 0 0 0 1 1 0
1 0.19 0.0361 0.45 1 0 0 1 2 1 0 0 1 0 0 0 2 1 1
0 0.19 0.0361 0.9 0 0 0 3 0 0 0 0 1 0 1 4 2 1 1
0 0.19 0.0361 0.15 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0
0 0.19 0.0361 0.45 0 0 0 2 5 1 1 0 1 0 0 0 3 1 2
1 0.19 0.0361 0.35 0 0 0 5 1 9 1 0 1 0 0 0 1 1 0
1 0.19 0.0361 0.55 0 0 0 4 0 2 0 0 1 0 0 0 0 0 0
1 0.19 0.0361 0.15 0 0 0 3 0 6 0 0 1 0 0 0 1 1 0
1 0.19 0.0361 0.65 1 0 0 2 0 5 0 0 1 0 0 0 1 0 1
0 0.19 0.0361 0.15 1 0 0 1 0 0 0 0 1 0 0 0 1 1 0
0 0.19 0.0361 0.45 0 0 0 1 0 0 0 0 1 0 0 0 1 1 0

   ... more lines ...   

The variable Doctorco represents doctor visits. Additional variables in the data set that you want to evaluate as determinants of doctor visits include Sex (coded as 0=male, 1=female), Age (age in years divided by 100), Illness (number of illnesses during the two-week interval, with five or more coded as five), Income (annual income in Australian dollars divided by 1,000), and Hscore (a score on a general health questionnaire, in which a high score indicates bad health).

You can load the DocVisit data set into your CAS session by naming your CAS engine libref in the first statement of the following DATA step:

data mycas.docvisit;
   set docvisit;
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 Poisson regression model in the Bayesian framework by using the No-U-Turn Sampler (NUTS) with multiple Markov chains. NUTS is specified in the SAMPLER= option in the BAYES statement, and you can set several values that control the sampler as arguments to the SAMPLER=NUTS option. The NCHAIN= option controls the number of Markov chains, and the NSAMPLE= option controls the size of the posterior sample per chain. PROC CNTSELECT automatically aggregates the results of each chain 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 cntselect data = mycas.docvisit;
   class sex illness;
   model doctorco=sex illness income hscore;
   bayes seed = 5317 nsample = 2000 nchain = 4
         sampler = nuts(ntune = 2000)
         statistics(bychain) = summary
         diagnostics(bychain) = (summary heidel);
run;

Output 10.3.1 shows the Bayesian estimation results. The "Posterior Summaries" table displays the estimates of the posterior mean and posterior quantiles that are computed using all four chains. The other posterior summary statistics tables that the STATISTICS= option produces would similarly be computed using all four chains. The "MCMC Diagnostic Summaries" table displays basic convergence diagnostics that are 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 that are performed for each chain. Specifically, for each parameter it shows the number of chains that passed each test, 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 would both be summarized similarly.

Detailed results for each chain are available in Output 10.3.2 as a result of the BYCHAIN options. The "Posterior Summaries by Chain" table displays the estimates of the posterior mean and posterior quantiles that are computed separately for each chain. Similarly, the "MCMC Diagnostic Summaries by Chain" table displays basic convergence diagnostics that are computed for each chain, and the "Heidelberger-Welch Diagnostics by Chain" table displays the Heidelberger-Welch diagnostics that are computed for each chain. For more information about the differences between the tables in Output 10.3.1 and Output 10.3.2, see the section Using Multiple Chains in Chapter 2, Introduction to Bayesian Analysis Procedures.