The coda Package
Contents
1. Martyn Plummer See Also meme meme List 32 time mcmc thin Thinning interval Description thin returns the interval between successive values of a time series thin x is equivalent to 1 frequency x This is a generic function Methods have been implemented for mcmc objects Usage thin x Arguments x Author s Martyn Plummer See Also time a regular time series a list of arguments time mcemc Time attributes for mcmc objects Description These are methods for memc objects for the generic time series functions Usage time x start x S3 me end x S3 me thin x Arguments x S3 met S3 met hod hod hod hod for class for class for class for class an meme or memc list object extra arguments for future methods cma cme cme emc traceplot 33 See Also time start frequency thin traceplot Trace plot of MCMC output Description Displays a plot of iterations vs sampled values for each variable in the chain with a separate plot per variable Usage traceplot x smooth TRUE col type ylab Arguments x An meme or memc list object smooth draw smooth line through trace plot col graphical parameter see par type graphical parameter see par ylab graphical parameter see par further graphical parameters Note You can call this functi2. TRUE bwf auto layout TRUE ask dev interactive 24 raftery diag Arguments x an object of class memc or meme list trace Plot trace of each variable density Plot density estimate of each variable smooth Draw a smooth line through trace plots bwf Bandwidth function for density plots auto layout Automatically generate output format ask Prompt user before each page of plots Further arguments Author s Martyn Plummer See Also densplot traceplot raftery diag Raftery and Lewis s diagnostic Description raftery diag is a run length control diagnostic based on a criterion of accuracy of estimation of the quantile q It is intended for use on a short pilot run of a Markov chain The number of iterations required to estimate the quantile g to within an accuracy of r with probability p is calculated Separate calculations are performed for each variable within each chain If the number of iterations in data is too small an error message is printed indicating the minimum length of pilot run The minimum length is the required sample size for a chain with no correlation between consecutive samples Positive autocorrelation will increase the required sample size above this minimum value An estimate I the dependence factor of the extent to which autocorrela tion inflates the required sample size is also provided Values of I larger than 5 indicate strong autocorrelation which may be due to a poor choice of start
3. crosscorr x Arguments x an meme or memc list object FALSI E 10 cumuplot Value A matrix or 3 d array containing the correlations See Also ecrosscorr plot au tocorr crosscorr plot Plot image of correlation matrix Description crosscorr plot provides an image of the correlation matrix for x If x is an mcmc list object then all chains are combined The range 1 1 is divided into a number of equal length categories given by the length of col and assigned the corresponding color By default topographic colours are used as this makes it easier to distinguish positive and negative correlations Usage crosscorr plot x col topo colors 10 Arguments x an meme or mcmc list object col color palette to use graphical parameters See Also erosscorr image topo colors cumuplot Cumulative quantile plot Description Plots the evolution of the sample quantiles as a function of the number of iterations Usage cumuplot x probs c 0 025 0 5 0 975 ylab lty c 2 1 lwd c 1 2 type 1 ask dev interactive auto layout TRUE col 1 densplot 11 Arguments x an memc object probs vector of desired quantiles ylab lty lwd type col graphical parameters auto layout If TRUE then set up own layout for plots otherwise use existing one ask If TRUE then prompt user before displaying each page of plots further graphical parameters Autho
4. Author s Deepayan Sarkar Deepayan Sarkar R project org varnames 37 See Also Lattice for a brief introduction to lattice displays and links to further documentation Examples data line xyplot line xyplot line 1 1 start 10 densityplot line start qqmath line start 10 levelplot line 2 acfplot line outer TRUE ll m o varnames Named dimensions of MCMC objects Description varnames returns the variable names and chanames returns the chain names or NULL if these are not set If allow null FALSE then NULL values will be replaced with canonical names Usage varnames x allow null TRUE chanames x allow null TRUE varnames x lt value chanames x lt value Arguments x an meme or mcmc list object allow null Logical argument that determines whether the function may return NULL value A character vector or NULL Value A character vector or NULL See Also meme meme list 38 window mcmc window mcmc Time windows for mcmc objects Description window mcmc is a method for mcmc objects which is normally called by the generic function window In addition to the generic parameters st art and end the additional parameter thin may be used to thin out the Markov chain Setting thin k selects every kth iteration starting with the first Note that the value of thin is absolute not relative The value supplied given
5. goodness of fit criteria based on stochastic processes Ann Math Statist 23 192 212 1952 Csorgo S and Faraway JJ The exact and asymptotic distributions of the Cramer von Mises statistic J Roy Stat Soc B 58 221 234 1996 HPDinterval 3 HPDinterval Highest Posterior Density intervals Description Create Highest Posterior Density HPD intervals for the parameters in an MCMC sample Usage HPDinterval obj prob 0 95 S3 method for class mcmc HPDinterval obj prob 0 95 S3 method for class mcmc list HPDinterval obj prob 0 95 Arguments obj The object containing the MCMC sample usually of class mcme or memc list prob A numeric scalar in the interval 0 1 giving the target probability content of the intervals The nominal probability content of the intervals is the multiple of 1 nrow obj nearest to prob Optional additional arguments for methods None are used at present Details For each parameter the interval is constructed from the empirical cdf of the sample as the shortest interval for which the difference in the ecdf values of the endpoints is the nominal probability Assuming that the distribution is not severely multimodal this is the HPD interval Value For an mcmc object a matrix with columns lower and upper and rows corresponding to the parameters The attribute Probability is the nominal probability content of the intervals A
6. indicating whether variables in x should be transformed to im prove the normality of the distribution If set to TRUE a log transform or logit transform as appropriate will be applied autoburnin a logical flag indicating whether only the second half of the series should be used in the computation If set to TRUE default and start x is less than end x 2 then start of series will be adjusted so that only second half of series is used Theory Gelman and Rubin 1992 propose a general approach to monitoring convergence of MCMC output in which two or more parallel chains are run with starting values that are overdispersed relative to the posterior distribution Convergence is diagnosed when the chains have forgotten their initial values and the output from all chains is indistinguishable The gelman diag diagnostic is ap plied to a single variable from the chain It is based a comparison of within chain and between chain variances and is similar to a classical analysis of variance There are two ways to estimate the variance of the stationary distribution the mean of the empirical variance within each chain W and the empirical variance from all chains combined which can be expressed as Bs 1 B W 7 mn 1B PIECES n n where B is the empirical between chain variance If the chains have converged then both estimates are unbiased Otherwise the first method will underestimate the variance since the individual chains have n
7. list of such matrices is returned for an mcmc list object Author s Douglas Bates Examples data line HPDinterval line 4 autocorr as ts mcmc Coerce mcmc object to time series Description as ts mcemc will coerce an memc object to a time series Usage as ts mcmc X Arguments x an memc object unused arguments for compatibility with generic as ts Author s Martyn Plummer See Also as ES autocorr Autocorrelation function for Markov chains Description autocorr calculates the autocorrelation function for the Markov chain memc obj at the lags given by lags The lag values are taken to be relative to the thinning interval if relat ive TRUE High autocorrelations within chains indicate slow mixing and usually slow convergence It may be useful to thin out a chain with high autocorrelations before calculating summary statistics a thinned chain may contain most of the information but take up less space in memory Re running the MCMC sampler with a different parameterization may help to reduce autocorrelation Usage E autocorr x lags c 0 1 5 10 50 relative TRUI Arguments x an memc object lags a vector of lags at which to calculate the autocorrelation relative a logical flag TRUE if lags are relative to the thinning interval of the chain or FALSE if they are absolute difference in iteration numbers autocorr diag 5 Value A vector or array containing th
8. mcmc list Replicated Markov Chain Monte Carlo Objects Description The function mcmc list is used to represent parallel runs of the same chain with different starting values and random seeds The list must be balanced each chain in the list must have the same iterations and the same variables Diagnostic functions which act on memc objects may also be applied to memc list objects In general the chains will be combined if this makes sense otherwise the diagnostic function will be applied separately to each chain in the list Since all the chains in the list have the same iterations a single time dimension can be ascribed to the list Hence there are time series methods time window start end frequency and thin for memc list objects Anmemc 1ist can be indexed as if it were a single memc object using the operator see exam ples below The operator selects a single mcmc object from the list Usage memc list as mcmc list x is mcemc list x Arguments a list of memc objects x an object that may be coerced to memc list Author s Martyn Plummer See Also meme Examples data line xl lt line 1 Select first chain x2 lt line 1 drop FALSE Select first var from all chains varnames x2 varnames line 1 TRUE memc subset 21 mcmc subset Extract or replace parts of MCMC objects Description These are methods for subsetting mcmc objects You can sele
9. message message displayed before prompting for user input what the type of what gives the type of data to be read lower lower limit of input for numeric input only upper upper limit of input for numeric input only answer in the input must correspond to one of the elements of the vector answer in if supplied default value assumed if user enters a blank line Value The value of the valid input When the default argument is specified a blank line is accepted as valid input and in this case read and check returns the value of default Note Since the function does not return a value until it receives valid input it extensively checks the conditions for consistency before prompting the user for input Inconsistent conditions will cause an error Author s Martyn Plummer read coda Read output files in CODA format Description read coda reads Markov Chain Monte Carlo output in the CODA format produced by Open BUGS and JAGS By default all of the data in the file is read but the arguments start end and thin may be used to read a subset of the data If the arguments given to start end or thin are incompatible with the data they are ignored Usage read coda output file index file start end thin quiet FALSE read jags file jags out start end thin quiet FALSE Arguments output file The name of the file containing the monitored output index file The name of the file containing the index sho
10. plots otherwise use existing one ask If TRUE then prompt user before displaying each page of plots graphical parameters See Also autocorr batchSE Batch Standard Error Description Effective standard deviation of population to produce the correct standard errors Usage batchSE x batchSize 100 Arguments x An meme or memc list object batchSize Number of observations to include in each batch Details Because of the autocorrelation the usual method of taking var x n overstates the precision of the estimate This method works around the problem by looking at the means of batches of the parameter If the batch size is large enough the batch means should be approximately uncorrelated and the normal formula for computing the standard error should work The batch standard error procedure is usually thought to be not as accurate as the time series meth ods used in summary and effectiveSize It is included here for completeness bugs2jags 7 Value A vector giving the standard error for each column of x Author s Russell Almond References Roberts GO 1996 Markov chain concepts related to sampling algorithms in Gilks WR Richard son S and Spiegelhalter DJ Markov Chain Monte Carlo in Practice Chapman and Hall 45 58 See Also spectrum0 ar effectiveSize summary mcmc bugs2jags Convert WinBUGS data file to JAGS data file Description bugs2jags converts a Wi
11. stationarity test and indicates that a longer MCMC run is needed If the stationarity test is passed the number of iterations to keep and the number to discard are reported The half width test calculates a 95 confidence interval for the mean using the portion of the chain which passed the stationarity test Half the width of this interval is compared with the estimate of the mean If the ratio between the half width and the mean is lower than eps the halfwidth test is passed Otherwise the length of the sample is deemed not long enough to estimate the mean with sufficient accuracy Theory The heidel diag diagnostic is based on the work of Heidelberger and Welch 1983 who com bined their earlier work on simulation run length control Heidelberger and Welch 1981 with the work of Schruben 1982 on detecting initial transients using Brownian bridge theory Note If the half width test fails then the run should be extended In order to avoid problems caused by sequential testing the test should not be repeated too frequently Heidelberger and Welch 1981 suggest increasing the run length by a factor I gt 1 5 each time so that estimate has the same reasonably large proportion of new data 18 mcmc References Heidelberger P and Welch PD A spectral method for confidence interval generation and run length control in simulations Comm ACM 24 233 245 1981 Heidelberger P and Welch PD Simulation run length control in the pres
12. subset 20 acf 4 35 acfplot trellisplots 32 acfplot mcemc list trellisplots 32 as array 19 as array mcmc list memc convert 18 as matrix 19 as matrix mcmc memc convert 18 as memc 19 as memc memc 17 as mecmc list memc list 19 as mecmc mcmc list memc convert 18 as ts 3 as ts mcmce 3 autocorr 3 4 5 9 autocorr diag 4 autocorr diag mcmc list autocorr diag 4 autocorr plot 4 5 batchSE 5 bugs2 Jags 6 40 chanames varnames 36 chanames lt varnames 36 coda options 7 codamenu 8 26 Cramer 1 crosscorr 8 9 crosscorr plot 9 9 cumuplot 9 density 10 densityplot mcmc trellisplots 32 densplot 10 23 32 display coda options coda options 7 dput 7 dump 7 ffectiveSize 6 11 end mcmc time mcmc 31 frequency 32 frequency mcmc time mcmc 31 gelman diag 11 74 gelman plot 13 3 geweke diag 14 15 geweke plot 15 5 glm 29 30 heidel diag 16 HPDinterval 2 HPDinterval mcmc list HPDinterval 2 image 9 is mcmc memc 17 is mcemc list mcmc list 19 Lattice 36 levelplot mcmc trellisplots 32 line 17 memc 17 19 22 26 30 36 meme convert 18 memc 1ist 18 19 21 22 26 30 36 memc subset 20 memcUpgrade 18 20 mcpar 21 menu 22 multi menu 21 nchain 22 niter nchain 22 INDEX nvar nchain 22 options 8 panel densityplot 35 panel xyplot 35 pcramer Cr
13. BK Ge S AOE REN N 18 MCMC CONVEIl 2 4 4 Wa Sees e ee Gee tre A AAA 19 Meme list cerrar ek hae SO Agony Botte od do ee Sore Sh tds 20 2 Cramer MGMGSUDSEE sb doch by Ste a de ed 21 memcUperade ss pesos eoa se RR Re eo 21 MCPAl arras A See S OMS SSS a AAA ENG 22 Mueni a SE Sear ers Sake amp be ew dt bees Coe Bako o 22 NERA 6 44 2 4 A a bd E E a aT 23 P O ICM pira a a 23 TAME Clas ci a ae A He a oy eS A a e a 24 co AA a ee ek es 25 A ey dered eo es MS Ae Ae etn ue ete Be ee we eae Bae eee 26 Tead codaunteractive cocos Hak ee REESE ORS a ERE HSS 27 read Openbues pco ope ps b rs bee G gue eb de Gee f 28 rejectionRate osc eo ca eR ER ER e eR Ra 29 SpecttumO co 4s a Dae RR a pea ee ee a eae ae ee 29 Spectra at ses eA a he OR ER A A Ra Oe a 30 summary MCMC 6 2g ea Rhea et ee ee ea eda YS 31 E eee Gc os AEA en ee eee OE Re ee ee 32 HMEMCMC essa e Ae we EER Se we HA Be Ss 32 traceplots 4 224 bb G44 we a Ge ba eR Dee ee ede 33 trellisplotS s mii ce i op a e a RR RE A Ae a Bae ae 33 VAINAMES aoe ee eo Be a Oe we Ea ee ee a eB a 37 WINdOWANCMG 2 00 a a RRE ae ee a ie ee ee ee ee 2 38 Index 39 Cramer The Cramer von Mises Distribution Description Distribution function of the Cramer von Mises distribution Usage pcramer q eps Arguments q vector of quantiles eps accuracy required Value pcramer gives the distribution function References Anderson TW and Darling DA Asymptotic theory of certain
14. The coda Package October 11 2006 Version 0 10 7 Date 2006 08 29 Title Output analysis and diagnostics for MCMC Author Martyn Plummer Nicky Best Kate Cowles Karen Vines Maintainer Martyn Plummer lt plummerGiarc fr gt Depends R gt 2 2 0 lattice Description Output analysis and diagnostics for Markov Chain Monte Carlo simulations License GPL Version 2 or later R topics documented Cramer aera BAe Oe eee RORY ER SSS BAe SY aE Oe BA EE a 8S 2 HPDintetvall ses p s e a ee Sodom ee ee a ON bee me Be a ee Se ee e e 3 as iS MEM 2454 24 ewe aa eed ehh 42 Dee SP ad S 4 AULOCOM gx sje e oe we EE won bee Beg dh le Rh we eS Bw we he wl 4 A LOCOTEdIAS 24 aa eee a Swe A a we ba ae 3 autoco plot ici Ba RG Be SAG Be eo Ge Raa 2 6 batchSE 6 25545 ce aa Gawd bee ee MA REE ORS Ee hoa L 6 DUESLIALS ceci ha ee Ate A um Bi ee o e ae tig d A es a dies 7 coda options ensk e 4 Bon Soe eM A aoe Sebi eo 6 Ee Se ae SR Soe 8 CO AMEN e 44 32 dee a BEd ba Se BS Se DEG A a 9 CIOSSCOT oos ses 4 5 eR ee ea oE a E a e 9 crosscom plot dis a ads 4 48 S804 A 10 CUMUPIOt 2 as e RE E A A SA A Se Sok es 10 densplot sca rare ee eb ea eee eS ee as 11 CHECHVESIZE cost 5 a We G Soe Bava ha 12 pelman di g evo ee a ee A a eh ae ee 4 12 gelman plot s pss ap A aa A a a a o as 14 PCWeke dag c aa oaae a i ER e a E ee ee Ee ale 15 geWwekKe plOt sir e e a aa a e e a a 16 Heidel Mae a ds e dd 17 line carrasco a a dl De AA E A 18 MEME gecs oS RE ERE EY A REESE
15. amer plot mcmc 10 18 22 32 plot mcmc list memc list 19 print 35 print mcmc memc 17 qqmath mcmc trellisplots 32 raftery diag 23 read and check 24 read coda 25 26 27 read coda interactive 26 26 read jags read coda 25 read openbugs 26 27 rejectionRate 28 rejectionRate mcmc list rejectionRate 28 spectrum 29 30 spectrum 28 30 spectrum0 ar 6 11 29 29 start 32 start mcmc time mcmc 31 summary mcmc 6 18 30 thin 18 31 32 37 thin memc time mcmc 31 time 31 32 time mcmc 31 topo colors 9 traceplot 23 32 trellisplots 32 ts 2 update 35 varnames 36 varnames lt varnames 36 window 37 window mecmc 18 20 37 xyplot 35 xyplot memc trellisplots 32
16. an meme or memc list object ignored present for consistency with generic for the mcmc list methods a logical flag to control whether multiple runs of a series are displayed in the same panel they are if FALSE not if TRUI 2 If 36 trellisplots specified in the mcmc methods this argument is ignored with a warning groups for the memc list methods a logical flag to control whether the underlying lattice call will be supplied a groups arguments indicating which run a data point originated from The panel function is responsible for handling such an argument and will usually differentiate runs within a panel by using different graphical parameters When outer FALSE the default of groups is TRUE if the corresponding default panel function is able to make use of such information When outer FALSE groups TRUE will be ignored with a warning aspect controls the physical aspect ratio of the panel See xyplot for details The default for these methods is chosen carefully check what the default plot looks like before changing this parameter default scales this parameter provides a reasonable default value of the scales parameter for the method It is unlikely that a user will wish to change this parameter Pass a value for scales see xyplot instead which will override values specified here type a Character vector that determines if lines points etc are drawn on the panel The default value
17. ardless of whether the variables were drawn separately or in a block In the latter case the values returned should be the same Usage rejectionRate x Arguments xX An memc or memc list object Details For the purposes of this function a rejection has occurred if the value of the time series is the same at two successive time points This test is done naively using and may produce problems due to rounding error Value A vector containing the rejection rates one for each variable Author s Russell Almond spectrum0 Estimate spectral density at zero Description The spectral density at frequency zero is estimated by fitting a glm to the low frequency end of the periodogram spect rum0 x length x estimates the variance of mean x Usage spectrum0 x max freq 0 5 order 1 max length 200 Arguments x A time series max freq The glm is fitted on the frequency range 0 max freq order Order of the polynomial to fit to the periodogram max length The data x is aggregated if necessary by taking batch means so that the length of the series is less than max length If this is set to NULL no aggregation occurs 30 spectrum0 ar Details The raw periodogram is calculated for the series x and a generalized linear model with family Gamma and log link is fitted to the periodogram The linear predictor is a polynomial in terms of the frequency The degree of the polynomial is determined by t
18. attr x title ylab or subset S3 method for class mcmc list qqmath x data outer FALSE groups outer aspect xy default scales list y list relation free prepanel prepanel qqmathline start 1 thin 1 main attr x title ylab Hi un or subset S3 method for class memc trellisplots xyplot x xyplot x 35 data outer layout c 1 ncol x default scales list y list relation free type l start 1 thin 1 ylab xlab Iteration number main attr x title or subset S3 method for class mcmc list data outer FALSE groups outer aspect xy layout c 1 ncol x 1 default scales list y list relation free type l start 1 thin 1 main attr x title ylab oe 4 subset data acfplot x S3 method for class mcmc acfplot x data outer prepanel panel type h aspect xy start 1 thin 1 lag max NULL ylab Autocorrelation xlab Lag main ater x Erte tyy ae d subset HH S3 method for class mcmc list acfplot x Arguments x data outer data outer FALSE groups outer prepanel panel type if groups b else h aspect xy start 1 thin 1 lag max NULL ylab Autocorrelation xlab Lag attr sp MEitliemy main E subset
19. ct iterations using the first dimension and variables using the second dimension Selecting iterations will return a vector or matrix not an meme object If you want to do row subsetting of an memc object and preserve its dimensions use the window function Subsetting applied to an memc 1ist object will simultaneously affect all the parallel chains in the object Usage x i j Arguments x An meme object i Row to extract j Column to extract See Also window mcmc mcmcUpgrade Upgrade mcmc objects in obsolete format Description In previous releases of CODA an mcmc object could be a single or multiple chains A new class mcmc list has now been introduced to deal with multiple chains and memc objects can only have data from a single chain Objects stored in the old format are now obsolete and must be upgraded Usage mcmcUpgrade x Arguments x an obsolete memc object Author s Martyn Plummer See Also meme 22 multi menu mcpar Mcpar attribute of MCMC objects Description The mcpar attribute of an MCMC object gives the start iteration the end iteration and the thinning interval of the chain It resembles the tsp attribute of time series t s objects Usage mcpar x Arguments x An mcmcm or memc list object See Also ts memc meme List multi menu Choose multiple options from a menu Description multi menu presents the user with a menu of choices
20. e autocorrelations Author s Martyn Plummer See Also act autocorr plot autocorr diag Autocorrelation function for Markov chains Description autocorr diag calculates the autocorrelation function for the Markov chain memc obj at the lags given by lags The lag values are taken to be relative to the thinning interval if relat ive TRUI Unlike autocorr if memc obj has many parmeters it only computes the autocorrelations with itself and not the cross correlations In cases where autocorr would return a matrix this function returns the diagonal of the matrix Hence it is more useful for chains with many parameters but may not be as helpful at spotting parameters E If memc obj is of class memc 11 st then the returned vector is the average autocorrelation across all chains Usage autocorr diag mcmc obj Arguments mame obj an object of class memc ormemc list optional arguments to be passed to autocorr Value A vector containing the autocorrelations Author s Russell Almond See Also autocorr ac autocorr plot 6 batchSE autocorr plot Plot autocorrelations for Markov Chains Description Plots the autocorrelation function for each variable in each chain in x Usage autocorr plot x lag max auto layout TRUE ask dev interactive Arguments x A Markov Chain lag max Maximum value at which to calculate acf auto layout If TRUE then set up own layout for
21. ence of an initial transient Opns Res 31 1109 44 1983 Schruben LW Detecting initialization bias in simulation experiments Opns Res 30 569 590 1982 line Simple linear regression example Description Sample MCMC output from a simple linear regression model given in the BUGS manual Usage data line Format An meme object Source Spiegelhalter D J Thomas A Best N G and Gilks W R 1995 BUGS Bayesian inference using Gibbs Sampling Version 0 5 MRC Biostatistics Unit Cambridge memc Markov Chain Monte Carlo Objects Description The function memc is used to create a Markov Chain Monte Carlo object The data are taken to be a vector or a matrix with one column per variable An memc object may be summarized by the summary function and visualized with the plot function MCMC objects resemble time series t s objects and have methods for the generic functions t ime start end frequency and window Usage mcmc data NA start 1 nd numeric 0 thin 1 as mcmc x is mcmc x memc convert Arguments data start end thin x Note 19 a vector or matrix of MCMC output the iteration number of the first observation the iteration number of the last observation the thinning interval between consecutive observations An object that may be coerced to an memc object The format of the memc class has changed between coda version 0 3 and 0 4 Older memc ob
22. etty printed using the display coda options function which groups the options into sections Available options are bandwidth Bandwidth function used when smoothing samples to produce density estimates De faults to Silverman s Rule of thumb combine corr Logical option that determines whether to combine multiple chains when calculating cross correlations combine plots Logical option that determines whether to combine multiple chains when plotting combine plots Logical option that determines whether to combine multiple chains when calculat ing summary statistics data saved For internal use only densplot Logical option that determines whether to plot a density plot when plot methods are called for memc objects digits Number of significant digits to use when printing fracl For Geweke diagnostic fraction to use from start of chain Defaults to 0 1 frac2 For Geweke diagnostic fraction to use from end of chain Default to 0 5 gr bin For Geweke Brooks plot number of iterations to use per bin gr max For Geweke Brooks plot maximum number of bins to use This option overrides gr bin halfwidth For Heidelberger and Welch diagnostic the target value for the ratio of half width to sample mean lowess Logical option that controls whether to plot a smooth line through a trace plot when plotting MCMC output q For Raftery and Lewis diagnostic the target quantile to be estimated r For Raftery and Lewis diagnostic t
23. he parameter order Value A list with the following values spec The predicted value of the spectral density at frequency zero Theory Heidelberger and Welch 1991 observed that the usual non parametric estimator of the spectral den sity obtained by smoothing the periodogram is not appropriate for frequency zero They proposed an alternative parametric method which consisted of fitting a linear model to the log periodogram of the batched time series Some technical problems with model fitting in their original proposal can be overcome by using a generalized linear model Batching of the data originally proposed in order to save space has the side effect of flattening the spectral density and making a polynomial fit more reasonable Fitting a polynomial of degree zero is equivalent to using the batched means method Note The definition of the spectral density used here differs from that used by spec pgram We con sider the frequency range to be between 0 and 0 5 not between 0 and frequency x 2 The model fitting may fail on chains with very high autocorrelation References Heidelberger P and Welch P D A spectral method for confidence interval generation and run length control in simulations Communications of the ACM Vol 24 pp233 245 1981 See Also spectrum spectrum0 ar glm spectrum0 ar Estimate spectral density at zero Description The spectral density at frequency zero is estimated by fitting a
24. he required precision s For Raftery and Lewis diagnostic the probability of obtaining an estimate in the interval q r qtr quantiles Vector of quantiles to print when calculating summary statistics for MCMC output trace Logical option that determines whether to plot a trace of the sampled output when plotting MCMC output user layout Logical option that determines whether current value of par mfrow should be used for plots TRUE or whether the optimal layout should be calculated FALSE codamenu 9 Usage coda options display coda options data FALSE stats FALSE plots FALSE diags Coda Options Coda Options Default Arguments data logical flag show data options stats logical flag show summary statistic options plots logical flag show plotting options diags logical flag show plotting options list of options See Also options codamenu Main menu driver for the coda package Description codamenu presents a simple menu based interface to the functions in the coda package It is designed for users who know nothing about the R language Usage codamenu Author s Kate Cowles Nicky Best Karen Vines Martyn Plummer crosscorr Cross correlations for MCMC output Description crosscorr calculates cross correlations between variables in Markov Chain Monte Carlo output If x is an mcmc list then all chains in x are combined before calculating the correlation Usage
25. ing value high posterior correlations or stickiness of the MCMC algorithm The number of burn in iterations to be discarded at the beginning of the chain is also calculated Usage raftery diag data q 0 025 r 0 005 s 0 95 converge eps 0 001 Arguments data an memc object q the quantile to be estimated r the desired margin of error of the estimate s the probability of obtaining an estimate in the interval q r q r converge eps Precision required for estimate of time to convergence read and check 25 Value A list with class raftery diag A print method is available for objects of this class the contents of the list are tspar The time series parameters of data params A vector containing the parameters r s and q Niters The number of iterations in data resmatrix A 3 d array containing the results M the length of burn in N the required sample size N min the minimum sample size based on zero autocorrelation and I M N Nmin the dependence factor Theory The estimated sample size for variable U is based on the process Z d U lt u where d is the indicator function and u is the qth quantile of U The process Z is derived from the Markov chain data by marginalization and truncation but is not itself a Markov chain However Z may behave as a Markov chain if it is sufficiently thinned out raftery diag calculates the smallest value of thinning interval k which makes the thinned chain Z behave as a Mar
26. ion Estimation of the effective sample size requires estimating the spectral density at frequency zero This is done by the function spectrum0 ar For a memc list object the effective sizes are summed across chains To get the size for each chain individually use lapply x effectiveSize Value A vector giving the effective sample size for each column of x See Also spectrum0 ar gelman diag Gelman and Rubin s convergence diagnostic Description The potential scale reduction factor is calculated for each variable in x together with upper and lower confidence limits Approximate convergence is diagnosed when the upper limit is close to 1 For multivariate chains a multivariate value is calculated that bounds above the potential scale reduction factor for any linear combination of the possibly transformed variables The confidence limits are based on the assumption that the stationary distribution of the variable under examination is normal Hence the transform parameter may be used to improve the normal approximation Usage gelman diag x confidence 0 95 transform FALSE autoburnin TRUI El gelman diag 13 Arguments x Anmemc list object with more than one chain and with starting values that are overdispersed with respect to the posterior distribution confidence the coverage probability of the confidence interval for the potential scale reduc tion factor transform a logical flag
27. jects will now cause is mcmc to fail with an appropriate warning message Obsolete memc objects can be upgraded with the mcmcUpgrade function Author s Martyn Plummer See Also mcemc list mcmcUpgrade thin window mcmc summary mcmc plot meme mcmc convert Conversions of MCMC objects Description These are methods for the generic functions as matrix as array andas memc as matrix strips the MCMC attributes from an mcmc object and returns a matrix If iters TRUE then a column is added with the iteration number For mcmc list objects the rows of multiple chains are concatenated and if chains TRUE a column is added with the chain number mcmc list objects can be coerced to 3 dimensional arrays with the as array function An mcemc list object with a single chain can be coerced to an memc object with as memc If the argument has multiple chains this causes an error Usage as matrix mcmc x iters as matrix mcmc list x ite as array mcmc list x drop Arguments x iters chains drop FALSE rs FALSE chains FALS sea r E E An memc or memc list object logical flag add column for iteration number logical flag add column for chain number if mcmc list logical flag if TRU E the result is coerced to the lowest possible dimension further arguments for future methods 20 memc list See Also as matrix as array as mcmc
28. kov chain The required sample size is calculated from this thinned sequence Since some data is thrown away the sample size estimates are conservative The criterion for the number of burn in iterations m to be discarded is that the conditional dis tribution of ZE given Zo should be within converge eps of the equilibrium distribution of the chain ZF Note raftery diag is based on the FORTRAN program gibbsit written by Steven Lewis and avail able from the Statlib archive References 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 Jn Practical Markov Chain Monte Carlo W R Gilks D J Spiegel halter and S Richardson eds London U K Chapman and Hall read and check Read data interactively and check that it satisfies conditions Description Input is read interactively and checked against conditions specified by the arguments what lower upper and answer in If the input does not satisfy all the conditions an appropriate error message is produced and the user is prompted to provide input This process is repeated until a valid input value is entered Usage read and check message what numeric lower upper answer in default 26 read coda Arguments
29. labelled from 1 to the number of choices The user may choose one or more options by entering a comma separated list A range of values wan may also be specified using the operator Mixed expressions such as 1 3 5 6 are permitted If allow zero is set to TRUE one can select 0 to exit without choosing an item Usage GI multi menu choices title header allow zero TRUI Arguments choices Character vector of labels for choices title Title printed before menu header Character vector of length 2 giving column titles allow zero Permit 0 as an acceptable response Value Numeric vector giving the numbers of the options selected or 0 if no selection is made Author s Martyn Plummer nchain 23 See Also menu nchain Dimensions of MCMC objects Description These functions give the dimensions of an MCMC object niter x returns the number of iterations nvar x returns the number of variables chain x returns the number of parallel chains Usage niter x nvar x nchain x Arguments x An meme or memc list object Value A numeric vector of length 1 See Also meme meme List plot mcmc Summary plots of memc objects Description plot mecmc summarizes an memc or memc list object with a trace of the sampled output and a density estimate for each variable in the chain Usage S3 method for class mcmc plot x trace TRUE density TRUE smooth
30. n frac2 fraction to use from end of chain nbins Number of segments pvalue p value used to plot confidence limits for the null hypothesis auto layout If TRUE then set up own layout for plots otherwise use existing one ask Prompt user before displaying each page of plots Graphical parameters Note The graphical implementation of Geweke s diagnostic was suggested by Steve Brooks See Also geweke diag heidel diag 17 heidel diag Heidelberger and Welch s convergence diagnostic Description heidel diag is a run length control diagnostic based on a criterion of relative accuracy for the estimate of the mean The default setting corresponds to a relative accuracy of two significant digits heidel diag also implements a convergence diagnostic and removes up to half the chain in order to ensure that the means are estimated from a chain that has converged Usage heidel diag x eps 0 1 pvalue 0 05 Arguments x eps Target value for ratio of halfwidth to sample mean pvalue significance level to use Details The convergence test uses the Cramer von Mises statistic to test 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 The latter outcome constitutes failure of the
31. n autoregressive model spect rum0 x length x estimates the variance of mean x Usage spectrum0 ar x summary mcmc 31 Arguments x A time series Details The ar function to fit an autoregressive model to the time series x For multivariate time series separate models are fitted for each column The value of the spectral density at zero is then given by a well known formula Value A list with the following values spec The predicted value of the spectral density at frequency zero order The order of the fitted model Note The definition of the spectral density used here differs from that used by spec pgram We con sider the frequency range to be between 0 and 0 5 not between 0 and frequency x 2 See Also spectrum spectrum0 glm summary mcmc Summary statistics for Markov Chain Monte Carlo chains Description summary mcmc produces two sets of summary statistics for each variable Mean standard deviation naive standard error of the mean ignoring autocorrelation of the chain and time series standard error based on an estimate of the spectral density at 0 Quantiles of the sample distribution using the quantiles argument Usage S3 method for class mcmc summary object quantiles c 0 025 0 25 0 5 0 75 0 975 Arguments object an object of class memc or memc list quantiles a vector of quantiles to evaluate for each variable a list of further arguments Author s
32. nBUGS data in the format called S Plus i e the format created by the dput function and writes it in dump format used by JAGS NB WinBUGS stores its arrays in row order This is different from R and JAGS which both store arrays in column order This difference is taken into account by bugs2 jags which will automati cally reorder the data in arrays without changing the dimension Usage bugs2jags infile outfile Arguments infile name of the input file outfile name of the output file Note If the input file is saved from WinBUGS it must be saved in plain text format The default format for files saved from WinBUGS is a binary compound document format with extension odc that cannot be read by bugs2jags Author s Martyn Plummer References Spiegelhalter DJ Thomas A Best NG and Lunn D 2003 WinBUGS version 1 4 user manual MRC Biostatistics Unit Cambridge UK 8 coda options See Also dput dump coda options Options settings for the codamenu driver Description coda options is a utility function that queries and sets options for the codamenu function These settings affect the behaviour of the functions in the coda library only when they are called via the codamenu interface The coda options function behaves just like the options function in the base library with the additional feature that coda options default TRUE will reset all options to the default values Options can be pr
33. ng Version 0 50 MRC Biostatistics Unit Cambridge See Also mcmc memc 1ist read coda codamenu 28 read openbugs read openbugs Read CODA output files produced by OpenBUGS Description read openbugs reads Markov Chain Monte Carlo output in the CODA format produced by OpenBUGS This is a convenience wrapper around the function read coda which allows you to read all the data output by OpenBUGS by specifying only the file stem Usage read openbugs stem start end thin quiet FALSE Arguments stem Character string giving the stem for the output files OpenBUGS produces files with names lt stem gt CODAindex txt lt stem gt CODAchain1 txt lt stem gt CODAchain2 txt start First iteration of chain end Last iteration of chain thin Thinning interval for chain quiet Logical flag If true a progress summary will be printed Value An object of class memc 1ist containing output from all chains Author s Martyn Plummer References Spiegelhalter DJ Thomas A Best NG and Lunn D 2004 WinBUGS User Manual Version 2 0 June 2004 MRC Biostatistics Unit Cambridge See Also read coda rejectionRate 29 rejectionRate Rejection Rate for Metropolis Hastings chains Description rejectionRate calculates the fraction of time that a Metropolis Hastings type chain rejected a proposed move The rejection rate is calculates separately for each variable in the memc obj argument irreg
34. on directly but it is more usually called by the plot mcmc function See Also densplot plot mcmc trellisplots Trellis plots for mcmc objects Description These methods use the Trellis framework as implemented in the lattice package to produce space conserving diagnostic plots from mcmc and memc list objects The xyplot meth ods produce trace plots The densityplot methods and qqmath methods produce empirical density and probability plots The levelplot method depicts the correlation of the series The acfplot methods plot the auto correlation in the series 34 trellisplots Usage S3 method for class mecmc densityplot x data outer aspect xy default scales list relation free start 1 thin 1 main attr x title xlab plot points rug r subset S3 method for class mcmc list densityplot x data outer FALSE groups outer aspect xy default scales list relation free start 1 thin 1 main attr x title xlab plot points rug let subset HH S3 method for class mecmc levelplot x data main attr x title start 1 thin 1 DAAT xlab ylab cuts 10 at col regions topo colors 100 subset S3 method for class mcmc qqmath x data outer aspect xy default scales list y list relation free prepanel prepanel qqmathline start 1 thin 1 main
35. ot had time to range all over the stationary distribution and the second method will overestimate the variance since the starting points were chosen to be overdispersed The convergence diagnostic is based on the assumption that the target distribution is normal A Bayesian credible interval can be constructed using a t distribution with mean jt Sample mean of all chains combined and variance A B V mn where m is the number of chains and degrees of freedom estimated by the method of moments o 2x V Var V Use of the t distribution accounts for the fact that the mean and variance of the posterior distribution are estimated The convergence diagnostic itself is Values substantially above 1 indicate lack of convergence If the chains have not converged Bayesian credible intervals based on the t distribution are too wide and have the potential to shrink by this factor if the MCMC run is continued 14 gelman plot Note The multivariate a version of Gelman and Rubin s diagnostic was proposed by Brooks and Gelman 1997 References Gelman A and Rubin DB 1992 Inference from iterative simulation using multiple sequences Statistical Science 7 457 511 Brooks SP and Gelman A 1997 General methods for monitoring convergence of iterative simu lations Journal of Computational and Graphical Statistics 7 434 455 See Also gelman plot gelman plot Gelman Rubin Brooks plot Description This
36. plot shows the evolution of Gelman and Rubin s shrink factor as the number of iterations increases Usage gelman plot x bin width 10 max bins 50 confidence 0 95 transform FALSE auto layout TRUE ask dev interactive col lty xlab ylab type Arguments x an memc object bin width Number of observations per segment excluding the first segment which always has at least 50 iterations max bins Maximum number of bins excluding the last one confidence Coverage probability of confidence interval transform Automatic variable transformation see gelman diag auto layout If TRUE then set up own layout for plots otherwise use existing one ask Prompt user before displaying each page of plots col graphical parameter see par lty graphical parameter see par xlab graphical parameter see par ylab graphical parameter see par type graphical parameter see par further graphical parameters geweke diag 15 Details The Markov chain is divided into bins according to the arguments bin width and max bins Then the Gelman Rubin shrink factor is repeatedly calculated The first shrink factor is calculated with observations 1 50 the second with observations 1 50 n where n is the bin width the third contains samples 1 50 2n and so on Theory A potential problem with gelman diag is that it may mis diagnose convergence if the shrink factor happens to be close to by chance By calculating
37. r s Arni Magnusson lt arnima u washington edu gt densplot Probability density function estimate from MCMC output Description Displays a plot of the density estimate for each variable in x calculated by the density function Usage densplot x show obs TRUE bwf main ylim Arguments x An memc or memc list object show obs Show observations along the x axis bwf Function for calculating the bandwidth If omitted the bandwidth is calculate by 1 06 times the minimum of the standard deviation and the interquartile range divided by 1 34 times the sample size to the negative one fifth power main Title See par ylim Limits on y axis See par Further graphical parameters Note You can call this function directly but it is more usually called by the plot mcmc function If a variable is bounded below at O or bounded in the interval 0 1 then the data are reflected at the boundary before being passed to the density function This allows correct estimation of a non zero density at the boundary See Also density plot mcmc 12 gelman diag ffectiveSiz Effective sample size for estimating the mean Description Sample size adjusted for autocorrelation Usage ffectiveSize x Arguments x An meme or memc list object Details For a time series x of length N the standard error of the mean is var x n where n is the effective sample size n Nonly when there is no autocorrelat
38. s for the methods are carefully chosen See panel xyplot for possible values thin an optional thinning interval that is applied before the plot is drawn start an optional value for the starting point within the series Values before the start ing point are considered part of the burn in of the series and dropped plot points Character argument giving the style in which points are added to the plot See panel densityplot for details layout a method specific default for the layout argument to the lattice functions xlab ylab main Used to provide default axis annotations and plot labels cuts at defines number and location of values where colors change col regions Color palette used lag max maximum lag for which autocorrelation is computed By default the value cho sen by acf is used prepanel panel suitable prepanel and panel functions for acfplot The prepanel function omits the lag 0 auto correlation which is always 1 from the range calculations other arguments passed to the lattice function Documentation of the corre sponding generics in the lattice package should be consulted for possible arguments subset indices of the subset of the series to plot The default is constructed from the start and thin arguments Value An object of class trellis The relevant update method can be used to update components of the object and the print method usually called by default will plot it on an appropriate plotting device
39. separate statistic is calculated for each variable in each chain 16 geweke plot References Geweke J Evaluating the accuracy of sampling based approaches to calculating posterior mo ments In Bayesian Statistics 4 ed JM Bernado JO Berger AP Dawid and AFM Smith Clarendon Press Oxford UK See Also geweke plot geweke plot Geweke Brooks plot Description If geweke diag indicates that the first and last part of a sample from a Markov chain are not drawn from the same distribution it may be useful to discard the first few iterations to see if the rest of the chain has converged This plot shows what happens to Geweke s Z score when succes sively larger numbers of iterations are discarded from the beginning of the chain To preserve the asymptotic conditions required for Geweke s diagnostic the plot never discards more than half the chain The first half of the Markov chain is divided into nbins 1 segments then Geweke s Z score is repeatedly calculated The first Z score is calculated with all iterations in the chain the second after discarding the first segment the third after discarding the first two segments and so on The last Z score is calculated using only the samples in the second half of the chain Usage geweke plot x fracl 0 1 frac2 0 5 nbins 20 pvalue 0 05 auto layout TRUE ask dev interactive Arguments x an memc object fracl fraction to use from beginning of chai
40. the shrink factor at several points in time gelman plot shows if the shrink factor has really converged or whether it is still fluctuating References Brooks S P and Gelman A 1998 General Methods for Monitoring Convergence of Iterative Simulations Journal of Computational and Graphical Statistics 7 p434 455 See Also gelman diag geweke diag Geweke s convergence diagnostic Description Geweke 1992 proposed a convergence diagnostic for Markov chains 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 are equal and Geweke s statistic has an asymptotically standard normal distribution The test statistic is a standard Z score the difference between the two sample means divided by its estimated standard error The standard error is estimated from the spectral density at zero and so takes into account any autocorrelation The Z score is calculated under the assumption that the two parts of the chain are asymptotically independent which requires that the sum of frac1 and frac2 be strictly less than 1 Usage geweke diag x fracl 0 1 frac2 0 5 Arguments x an memc object fracl fraction to use from beginning of chain frac2 fraction to use from end of chain Value Z scores for a test of equality of means between the first and last parts of the chain A
41. to the parameter thin must be a multiple of thin x Values of start end and thin which are inconsistent with x are ignored but a warning message is issued Usage HH S3 method for class mcmc window x start end thin Arguments x an memc object start the first iteration of interest end the last iteration of interest thin the required interval between successive samples futher arguments for future methods See Also window thin Index Topic array crosscorr 8 mcmc convert 18 Topic datasets line 17 Topic distribution Cramer 1 Topic file bugs2 Jags 6 read coda 25 read coda interactive 26 read openbugs 27 Topic hplot autocorr plot 5 crosscorr plot 9 cumuplot 9 densplot 10 gelman plot 13 geweke plot 15 plot mcmc 22 traceplot 32 trellisplots 32 Topic htest gelman diag 11 geweke diag 14 heidel diag 16 HPDinterval 2 raftery diag 23 Topic manip varnames 36 Topic multivariate crosscorr 8 Topic ts as ts mcemc 3 autocorr 3 autocorr diag 4 batchSE 5 ffectiveSize 11 meme 17 mcmc list 19 mcmc subset 20 mcmcUpgrade 20 mcpar 21 39 nchain 22 rejectionRate 28 spectrum 28 spectrum0 ar 29 thin 31 time mcmc 31 window mcmc 37 Topic univar HPDinterval 2 summary mcmc 30 Topic utilities as ts mcmc 3 coda options 7 codamenu 8 multi menu 21 read and check 24 Coda Options coda options 7 20 memc mcmc
42. wing which rows of the output file correspond to which variables file For JAGS output the name of the output file The extension out may be omitted There must be a corresponding ind file with the same file stem start First iteration of chain end Last iteration of chain thin Thinning interval for chain quiet Logical flag If true a progress summary will be printed read coda interactive 27 Value An object of class memc containing a representation of the data in the file Author s Karen Vines Martyn Plummer References Spiegelhalter DJ Thomas A Best NG and Gilks WR 1995 BUGS Bayesian inference Using Gibbs Sampling Version 0 50 MRC Biostatistics Unit Cambridge See Also mcmc read coda interactive read openbugs read coda interactive Read CODA output files interactively Description read coda interactive reads Markov Chain Monte Carlo output in the format produced by the classic BUGS program No arguments are required Instead the user is prompted for the required information Usage read coda interactive Value An object of class memc 1ist containing a representation of the data in one or more BUGS output files Note This function is normally called by the codamenu function but can also be used on a stand alone basis Author s Nicky Best Martyn Plummer References Spiegelhalter DJ Thomas A Best NG and Gilks WR 1995 BUGS Bayesian inference Using Gibbs Sampli
Download Pdf Manuals
Related Search