The S-Plus code for implementing raftery:lewis:92a's convergence diagnostic in CODA is based on the gibbsit function contributed to the Statlib archive by Steven Lewis (e-mail statlib@lib.stat.cmu.edu)
raftery:lewis:92a's method applies to single chains. It is intended both to detect convergence to the stationary distribution and to provide bounds for the accuracy of the estimated quantiles of functions of variables of interest. The user must specify the quantile to be estimated (the CODA default is the 2.5th percentile), the desired degree of accuracy for the estimate of this quantile (the CODA default is 
0.005) and the required probability of attaining this degree of accuracy (the CODA default is 0.95). The CODA output then reports Nmin -- the minimum number of iterations that would be needed to estimate the specified quantile to the desired precision if the samples in the chain were independent. This is a theoretical value based on the binomial variance and provides a lower bound for the run-length of the Gibbs sampler. (Note however that in the unusual event that consecutive samples in the output chain are negatively correlated, fewer than Nmin iterations will be needed -- for example, see the line output below). Nmin will increase as the required probability and degree of accuracy increase. Somewhat counter-intuitively, it also will be larger when estimating quantiles close to the median compared to more extreme quantiles. If the chain length of the BUGS output is less than Nmin for the quantile, accuracy and probability values currently specified in CODA, the program will not compute the Raftery & Lewis diagnostics, but will issue an error message stating the value of Nmin. If sufficient BUGS iterations are available, CODA will report N -- the total number of iterations that should be run for each variable; M -- the number of initial iterations to discard as the `burn-in' ; and k -- the thinning interval to be used. The final column in the CODA output reports I = 
. This measures the increase in number of iterations needed to reach convergence due to dependence between the samples in the chain. Values of I much greater than 1.0 indicate high within-chain correlations and probable convergence failure (raftery:lewis:92b suggest that I > 5.0 often indicates problems). In this case, reparameterization of the model is advised.
Selecting Raftery & Lewis (option 3:) from the CODA Diagnostics Menu produces the following output for the line example:
RAFTERY AND LEWIS CONVERGENCE DIAGNOSTIC: ========================================= Iterations used = 1:200 Thinning interval = 1 Sample size per chain = 200 Quantile = 0.025 Accuracy = +/- 0.005 Probability = 0.95 Chain: line1 ============ ******* Error: ******* Nmin = 3746; Available chain length is 200 Re-run BUGS program for at least 3746 iterations OR reduce accuracy, probability and/or quantile to be estimated Chain: line2 ============ ******* Error: ******* Nmin = 3746; Available chain length is 200 Re-run BUGS program for at least 3746 iterations OR reduce accuracy, probability and/or quantile to be estimated
This implies that in order to estimate the 2.5% quantile of each parameter to within 0.005 with 95% probability, a minimum of 3746 iterations of the Gibbs sampler are needed. However, if we relax the precision and degree of certainty required for the estimate of this quantile to 0.02 and 90% respectively (see §5.7.3 for details of how to change these parameters in CODA), we obtain the following diagnostics:
RAFTERY AND LEWIS CONVERGENCE DIAGNOSTIC: ========================================= Iterations used = 1:200 Thinning interval = 1 Sample size per chain = 200 Quantile = 0.025 Accuracy = +/- 0.02 Probability = 0.9 Chain: line1 ============ -+----------+----------------------------------------------------+- | | Thin Burn-in Total Lower bound Dependence | | VARIABLE | (k) (M) (N) (Nmin) factor (I) | | ======== | ==== ======= ===== =========== ========== | | | | | alpha | 1 5 244 165 1.48 | | beta | 1 2 131 165 0.794 | | sigma | 1 2 160 165 0.97 | | | | -+----------+----------------------------------------------------+- Chain: line2 ============ -+----------+----------------------------------------------------+- | | Thin Burn-in Total Lower bound Dependence | | VARIABLE | (k) (M) (N) (Nmin) factor (I) | | ======== | ==== ======= ===== =========== ========== | | | | | alpha | 1 2 160 165 0.97 | | beta | 1 2 160 165 0.97 | | sigma | 1 2 160 165 0.97 | | | | -+----------+----------------------------------------------------+-
The small burn-in values reported above suggest that both chains converge almost immediately for each monitored variable in the line example. For alpha in chain line1, the first 5 iterations should be discarded, and an estimated 44 additional BUGS iterations (on top of the 200 already obtained) performed in order to estimate the 2.5th percentile of the posterior distribution to the specified accuracy and probability. For the remaining variables and chains, 200 iterations are more than sufficient to estimate this quantile to the required precision, with only the first 2 iterations needing to be discarded. Note that N < Nmin for each of the latter variables. This is reflected by I < 1, and indicates negative correlations between consecutive iterates in the BUGS output.