convdiag_bmrm {bayesMRM}

Convergence Diagnostics on MCMC samples in bmrm

Description

Compute convergence diagnostics of Geweke (1992), Heidelberger and Welch (1983), Raftery and Lewis(1992).

Usage

convdiag_bmrm(x , var="P", convdiag="geweke",print=TRUE,...)

Arguments

x

an object of class bmrm, the output of the bmrm function

var

name of a variable to which convergence disagnostics apply. It should be one of "A" (source contribution matrix), "P" (source composition or profile matrix), "Sigma" (error variance).

convdiag

vector of convergence diagnostic methods. It should be any subvector of ("geweke", "heidel","raftery" ) (default="geweke").

print

TRUE/FALSE, print convergence diagnostics results (default=TRUE)

...

arguments to be passed to methods

Details

Geweke's convergence diagnostic for Markov chains is based on a test for equality of the means of the first and last part of a Markov chain (by default the first 10% and the last 50%). If the samples are drawn from the stationary distribution of the chain, the two means should be equal and Geweke's statistic has an asymptotically standard normal distribution. We use the function geweke.diag in coda package (with default option) which provides the test statistics (standard Z-scores) and the upper bound of and p-values.

Heidelberger and Welch's convergence diagnostic tests the null hypothesis that the sampled values come from a stationary distribution. The test is successively applied, firstly to the whole chain, then after discarding the first 10%, 20%, ... of the chain until either the null hypothesis is accepted, or 50% of the chain has been discarded. We use the function heidel.diag (with default option) which provides the staionary test results and p-values.

Raftery and Lewis's diagnostic estimates the minimum number of iterations, burn-in, thinning interval for zero autocorrelation, satisfying specified conditions regarding quantile q of parameters of interest. The conditions are specified by a posterior quantile q of parameters, an acceptable tolerance (accuracy) r for q, a probability s of being within the interval q-r, q+r. We use the function raftery.diag (with default option).

Value

A list of convergence diagnostics results

convdiag

selected convergence diagnostic methods

geweke

Geweke's z-scores and p-values if convdiag includes "geweke", NULL if convdiag does not include "geweke"

heidel

Heidelberger and Welch's stationary test results and p-values if convdiag includes "heidel"; NULL if convdiag does not include "heidel"

raftery

Raftery and Lewis's estimates of burn-in, minimum number of iterations, and thinning if convdiag includes "raftery"; NULL if convdiag does not include "raftery"

References

Geweke, J.(1992) Evaluating the accuracy of sampling-based approaches to calculating posterior moments. In Bayesian Statistics 4 (ed JM Bernado, JO Berger, AP Dawid and AFM Smith). Clarendon Press.

Heidelberger P, and Welch PD. (1981) A spectral method for confidence interval generation and run length control in simulations. Comm. ACM. 24, 233-245.

Heidelberger P. and Welch PD.(1983) Simulation run length control in the presence of an initial transient. Opns Res., 31, 1109-44,Oxford, UK.

Plummer, M., Best, N., Cowles, K. and Vines K. (2006) CODA: Convergence Diagnosis and Output Analysis for MCMC, R News, Vol 6, pp. 7-11.

Raftery, A.E. and Lewis, S.M. (1992). One long run with diagnostics: Implementation strategies for Markov chain Monte Carlo. Statistical Science, 7, 493-497.

Raftery, A.E. and Lewis, S.M. (1995). The number of iterations, convergence diagnostics and generic Metropolis algorithms. In Practical Markov Chain Monte Carlo (W.R. Gilks, D.J. Spiegelhalter and S. Richardson, eds.). London, U.K.: Chapman and Hall.

Examples

data(Elpaso)
Y=Elpaso$Y ; muP=Elpaso$muP
q=nrow(muP)
out.Elpaso <- bmrm(Y,q,muP, nAdapt=1000,nBurnIn=5000,nIter=5000,nThin=1)
conv1<-convdiag_bmrm(out.Elpaso,var="P",convdiag="raftery" )
conv2<-convdiag_bmrm(out.Elpaso,var="A", convdiag="geweke")
conv3<-convdiag_bmrm(out.Elpaso,var="Sigma", convdiag=c("geweke","heidel"))
conv4<-convdiag_bmrm(out.Elpaso,var="Sigma", convdiag=c("geweke","heidel", "raftery"))