Convergence in MCMC

The MI Procedure

Convergence in MCMC

The theoretical convergence of the MCMC process has been explored under various conditions, as described in Schafer (1997, p. 70). However, in practice, verification of convergence is not a simple matter and cannot be easily implemented in the MI procedure.

The parameters used in the imputation step for each iteration can be saved in an output data set with the OUTITER= option. These include the means, standard deviations, covariances, the worst linear function, and observed-data LR statistics. You can then monitor the convergence in a single chain by displaying time-series plots and autocorrelations for those parameter values (Schafer 1997, p. 120). The time-series and autocorrelation function plots for parameters such as variable means, covariances, and the worst linear function can be displayed by specifying the TIMEPLOT and ACFPLOT option.

You can apply EM to a bootstrap sample to obtain overdispersed starting values for multiple chains (Gelman and Rubin 1992). This provides a conservative estimate of the number of iterations needed before each imputation.

The next four subsections provide useful statistics and plots that can be used to check the convergence of the MCMC process.

LR Statistics

You can save the observed-data likelihood ratio (LR) statistic in each iteration with the LR option in the OUTITER= data set. The statistic is based on the observed-data likelihood with parameter values used in the iteration and the observed-data maximum likelihood derived from the EM algorithm.

In each iteration, the LR statistic is given by

$-2 \, {\rm log} \, ( \frac { f( \hat{{\theta}}_{i}) } { f( \hat {{\theta}}) } )$

where $f( \hat {{\theta}})$ is the observed-data maximum likelihood derived from the EM algorithm and $f( \hat{{\theta}}_{i})$ is the observed-data likelihood for $\hat{{\theta}}_{i}$ used in the iteration.

Similarly, you can also save the observed-data LR posterior mode statistic for each iteration with the LR_POST option. This statistic is based on the observed-data posterior density with parameter values used in each iteration and the observed-data posterior mode derived from the EM algorithm for posterior mode.

For large samples, these LR statistics tends to be approximately ${\chi}^2$ distributed with degrees of freedom equal to the dimension of ${{\theta}}$ (Schafer 1997, p. 131). For example, with a large number of iterations, if the values of the LR statistic do not behave like a random sample from the described ${\chi}^2$ distribution, then there is evidence that the MCMC process has not converged.

Worst Linear Function of Parameters

The worst linear function (WLF) of parameters (Schafer 1997, pp. 129-131) is a scalar function of parameters mu and ${\Sigma}$ that is "worst" in the sense that its function values converge most slowly among parameters in the MCMC process. The convergence of this function is evidence that other parameters are likely to converge as well.

For linear functions of parameters ${{\theta}}=( {mu}, {\Sigma})$ ,a worst linear function of ${{\theta}}$ has the highest asymptotic rate of missing information. The function can be derived from the iterative values of ${{\theta}}$ near the posterior mode in the EM algorithm. That is, an estimated worst linear function of ${{\theta}}$ is

$w( {{\theta}} )=v' \, ( {{\theta}} - \hat{{\theta}} )$

where $\hat{{\theta}}$ is the posterior mode and the coefficients $v=\hat{{\theta}}_{(-1)} - \hat{{\theta}}$ is the difference between the estimated value of ${{\theta}}$ one step prior to convergence and the converged value $\hat{{\theta}}$ .

You can display the coefficients of the worst linear function, v, by specifying the WLF option in the MCMC statement. You can save the function value from each iteration in an OUTITER= data set by specifying the WLF option in the OUTITER option. You can also display the worst linear function values from iterations in an autocorrelation plot or a time-series plot by specifying WLF as an ACFPLOT or TIMEPLOT option, respectively.

Note that when the observed-data posterior is nearly normal, the WLF is one of the slowest functions to approach stationarity. When the posterior is not close to normal, other functions may take much longer than the WLF to converge, as described in Schafer (1997, p.130).

Time-Series Plot

A time-series plot for a parameter ${\xi}$ is a scatter plot of successive parameter estimates ${\xi}_i$ against the iteration number i. The plot provides a simple way to examine the convergence behavior of the estimation algorithm for ${\xi}$ .Long-term trends in the plot indicate that successive iterations are highly correlated and that the series of iterations has not converged.

You can display time-series plots for the worst linear function, the variable means, variable variances, and covariances of variables. You can also request logarithmic transformations for positive parameters in the plots with the LOG option. When a parameter value is less than or equal to zero, the value is not displayed in the corresponding plot.

By default, the MI procedure uses the plus sign (+) as the plot symbol to display the points with a height of one (percentage screen unit) in a time-series plot. You can use the SYMBOL=, CSYMBOL=, and HSYMBOL= options to change the shape, color, and height of the plot symbol.

By default, the plot title "Time-Series Plot" is displayed in a time-series plot. You can request another title by using the TITLE= option in TIMEPLOT. When another title is also specified in a TITLE statement, this title is displayed as the main title and the plot title is displayed as a subtitle in the plot.

You can use options in the GOPTIONS statement to change the color and height of the title. Refer to the chapter "The SAS/GRAPH Statements" in SAS/GRAPH Software: Reference, Version 8 for a description of title options. See Example 9.6 for a usage of the time-series plot.

Autocorrelation Function Plot

To examine relationships of successive parameter estimates ${\xi}$ ,the autocorrelation function (ACF) can be used. For a stationary series, ${\xi}_i, i \ge 1$ , in time series data, the autocorrelation function at lag k is

${\rho}_{k}=\frac{ {\rm Cov}({\xi}_i, {\xi}_{i+k}) } { {\rm Var}({\xi}_i) }$

The sample k_th order autocorrelation is computed as

$r_{k}=\frac{\sum_{i=1}^{n-k} ({\xi}_i-{\overline {\xi}}) ({\xi}_{i+k}-{\overline {\xi}}) } {\sum_{i=1}^n ({\xi}_i-{\overline {\xi}})^2}$

You can display autocorrelation function plots for the worst linear function, the variable means, variable variances, and covariances of variables. You can also request logarithmic transformations for parameters in the plots with the LOG option. When a parameter has values less than or equal to zero, the corresponding plot is not created.

You specify the maximum number of lags of the series with the NLAG= option. The autocorrelations at each lag less than or equal to the specified lag are displayed in the graph. In addition, the plot also displays approximate 95% confidence limits for the autocorrelations. At lag k, the confidence limits indicate a set of approximate 95% critical values for testing the hypothesis ${\rho}_{j}=0, j \ge k.$

By default, the MI procedure uses the star sign (*) as the plot symbol to display the points with a height of one (percentage screen unit) in the plot, a solid line to display the reference line of zero autocorrelation, vertical line segments to connect autocorrelations to the reference line, and a pair of dashed lines to display approximately 95% confidence limits for the autocorrelations.

You can use the SYMBOL=, CSYMBOL=, and HSYMBOL= options to change the shape, color, and height of the plot symbol, and the CNEEDLES= and WNEEDLES= options to change the color and width of the needles. You can also use the LREF=, CREF=, and WREF= options to change the line type, color, and width of the reference line. Similarly, you can use the LCONF=, CCONF=, and WCONF= options to change the line type, color, and width of the confidence limits.

By default, the plot title "Autocorrelation Plot" is displayed in a autocorrelation function plot. You can request another title by using the TITLE= option in ACFPLOT. When another title is also specified in a TITLE statement, this title is displayed as the main title and the plot title is displayed as a subtitle in the plot.

Chapter Contents
Previous
Next
Top