MCMC Diagnostic Plots

Description

Markov chain Monte Carlo diagnostic plots. The purpose of the package is to combine existing tools from the coda and lattice packages, and make it easy to adjust graphical details.

Details

Diagnostic plots:

Posterior plots:

Examples:

Note

browseVignettes() shows a vignette with all the example plots.

The plot functions assume that MCMC results are stored either as a plain numeric vector (single chain) or in a data.frame (multiple chains). The mcmc class is also supported.

Author(s)

Arni Magnusson and Ian Stewart.

References

Fournier, D. A., Skaug, H. J., Ancheta, J., Ianelli, J., Magnusson, A., Maunder, M. N., Nielsen, A. and Sibert, J. (2012) AD Model Builder: using automatic differentiation for statistical inference of highly parameterized complex nonlinear models. Optimization Methods and Software, 27, 233–249.

Magnusson, A., Punt, A. E. and Hilborn, R. (2013) Measuring uncertainty in fisheries stock assessment: the delta method, bootstrap, and MCMC. Fish and Fisheries, 14, 325–342.

See Also

The coda package is a suite of diagnostic functions and plots for MCMC analysis, many of which are used in plotMCMC.

Many plotMCMC graphics are trellis plots, rendered with the lattice package.

The functions Args and ll (package gdata) can be useful for browsing unwieldy functions and objects.

Plot MCMC Autocorrelation

Description

Plot Markov chain Monte Carlo autocorrelation over a range of lag values. This is a diagnostic plot for deciding whether a chain needs further thinning.

Usage

plotAuto(mcmc, thin=1, log=FALSE, base=10, main=NULL, xlab="Lag",
         ylab="Autocorrelation", lty=1, lwd=1, col="black", ...)

Arguments

Value

Null, but a plot is drawn on the current graphics device.

Note

The Args function from the gdata package is recommended for reviewing the arguments, instead of args.

See Also

autocorr.plot is the underlying plotting function, and window.mcmc is used to optionally thin the chain(s).

plotTrace, plotAuto, plotCumu, and plotSplom are diagnostic plots.

plotDens and plotQuant are posterior plots.

plotMCMC-package gives an overview of the package.

Examples

plotAuto(xpar$R0)
plotAuto(xpar$R0, thin=10)
plotAuto(xpar, lag.max=50, ann=FALSE, axes=FALSE)

Plot MCMC Cumulative Quantiles

Description

Plot Markov chain Monte Carlo cumulative quantiles. This is a diagnostic plot for deciding whether the chain has converged.

Usage

plotCumu(mcmc, probs=c(0.025,0.975), div=1, log=FALSE, base=10,
         main=NULL, xlab="Iterations", ylab="Value", lty.median=1,
         lwd.median=2, col.median="black", lty.outer=2, lwd.outer=1,
         col.outer="black", ...)

Arguments

Value

Null, but a plot is drawn on the current graphics device.

Note

The Args function from the gdata package is recommended for reviewing the arguments, instead of args.

See Also

cumuplot is the underlying plotting function, and quantile is called iteratively to calculate the cumulative quantiles.

plotTrace, plotAuto, plotCumu, and plotSplom are diagnostic plots.

plotDens and plotQuant are posterior plots.

plotMCMC-package gives an overview of the package.

Examples

plotCumu(xpar$R0, main="R0")
plotCumu(xpar$cSfull, main="cSfull")
plotCumu(xpar, probs=c(0.25,0.75), ann=FALSE, axes=FALSE)

Plot MCMC Density

Description

Plot Markov chain Monte Carlo density. This is an approximation of the posterior probability density function.

Usage

plotDens(mcmc, probs=c(0.025,0.975), points=FALSE, axes=TRUE,
         same.limits=FALSE, between=list(x=axes,y=axes), div=1,
         log=FALSE, base=10, main=NULL, xlab=NULL, ylab=NULL,
         cex.main=1.2, cex.lab=1, cex.axis=0.8, cex.strip=0.8,
         col.strip="gray95", las=0, tck=0.5, tick.number=5,
         lty.density=1, lwd.density=3, col.density="black",
         lty.median=2, lwd.median=1, col.median="darkgray", lty.outer=3,
         lwd.outer=1, col.outer="darkgray", pch="|", cex.points=1,
         col.points="black", plot=TRUE, ...)

Arguments

Value

When plot=TRUE, a trellis plot is drawn and a data frame is returned, containing the data used for plotting. When plot=FALSE, a trellis object is returned.

Note

The Args function from the gdata package is recommended for reviewing the arguments, instead of args.

See Also

xyplot and panel.densityplot are the underlying drawing functions, and link[coda]{densplot} is a similar non-trellis plot.

plotTrace, plotAuto, plotCumu, and plotSplom are diagnostic plots.

plotDens and plotQuant are posterior plots.

plotMCMC-package gives an overview of the package.

Examples

plotDens(xbio$"2004", points=TRUE, div=1000, main="2004\n",
         xlab="Biomass age 4+ (kt)", tick.number=6, strip=FALSE)
plotDens(xpar, xlab="Parameter value", ylab="Posterior density\n")

Plot MCMC Quantiles

Description

Plot quantiles of multiple Markov chain Monte Carlo chains, using bars, boxes, or lines.

Usage

plotQuant(mcmc, style="boxes", probs=c(0.025,0.975), axes=TRUE,
          names=NULL, ylim=NULL, yaxs="i", div=1, log=FALSE, base=10,
          main=NULL, xlab=NULL, ylab=NULL, cex.axis=0.8, las=1,
          tck=-0.015, tick.number=8, lty.median=1*(style!="bars"),
          lwd.median=1+1*(style!="boxes"), col.median="black",
          lty.outer=1+2*(style=="lines"), lwd.outer=1,
          col.outer="black", pch=16, cex=0.8, col="black",
          boxfill="darkgray", boxwex=0.5, staplewex=0.5, sfrac=0.005,
          mai=c(0.8,1,1,0.6),
          mgp=list(bottom=c(2,0.4,0),left=c(3,0.6,0),top=c(0,0.6,0),
          right=c(0,0.6,0)), ...)

Arguments

Value

List containing:

Note

With style="boxes", the quartiles are shown as boxes.

The Args function from the gdata package is recommended for reviewing the arguments, instead of args.

See Also

bxp, plotCI, and matplot are the underlying drawing functions.

plotTrace, plotAuto, plotCumu, and plotSplom are diagnostic plots.

plotDens and plotQuant are posterior plots.

plotMCMC-package gives an overview of the package.

Examples

plotQuant(xrec, names=substring(names(xrec),3), div=1000, xlab="Year",
          ylab="Recruitment (million one-year-olds)")
plotQuant(xbio, div=1000, xlab="Year", ylab="Biomass age 4+ (kt)")
plotQuant(xbio, style="bars", div=1000, sfrac=0, xlab="Year",
          ylab="Biomass age 4+ (kt)")
plotQuant(xbio, style="lines", div=1000, xlab="Year",
          ylab="Biomass age 4+ (kt)")
plotQuant(xpro, axes=1:2, div=1000, xlab="Year",
          ylab="Biomass age 4+ (kt)")

Plot MCMC Scatterplot Matrix

Description

Plot scatterplots of multiple Markov chain Monte Carlo chains. This is a diagnostic plot for deciding whether parameters are confounded. When parameter estimates are highly dependent on each other, it may undermine conclusions based on MCMC results of that model.

Usage

plotSplom(mcmc, axes=FALSE, between=0, div=1, log=FALSE, base=10, ...)

Arguments

Value

Null, but a plot is drawn on the current graphics device.

Note

The Args function from the gdata package is recommended for reviewing the arguments, instead of args.

See Also

pairs is the underlying drawing function, and splom is a similar trellis plot.

plotTrace, plotAuto, plotCumu, and plotSplom are diagnostic plots.

plotDens and plotQuant are posterior plots.

plotMCMC-package gives an overview of the package.

Examples

plotSplom(xpar, pch=".")
plotSplom(xpro, axes=TRUE, between=1, div=1000, main="Future biomass",
          cex.labels=1.5, pch=".", cex=3)

Plot MCMC Traces

Description

Plot Markov chain Monte Carlo traces. This is a diagnostic plot for deciding whether a chain shows unwanted trends.

Usage

plotTrace(mcmc, axes=FALSE, same.limits=FALSE,
          between=list(x=axes,y=axes), div=1, span=1/4, log=FALSE,
          base=10, main=NULL, xlab=NULL, ylab=NULL, cex.main=1.2,
          cex.lab=1, cex.axis=0.8, cex.strip=0.8, col.strip="gray95",
          las=0, tck=0.5, tick.number=5, lty.trace=1, lwd.trace=1,
          col.trace="gray", lty.median=1, lwd.median=1,
          col.median="black", lty.loess=2, lwd.loess=1,
          col.loess="black", plot=TRUE, ...)

Arguments

Value

When plot=TRUE, a trellis plot is drawn and a data frame is returned, containing the data used for plotting. When plot=FALSE, a trellis object is returned.

Note

The Args function from the gdata package is recommended for reviewing the arguments, instead of args.

See Also

xyplot and panel.loess are the underlying drawing functions, and traceplot is a similar non-trellis plot.

plotTrace, plotAuto, plotCumu, and plotSplom are diagnostic plots.

plotDens and plotQuant are posterior plots.

plotMCMC-package gives an overview of the package.

Examples

plotTrace(xpar, xlab="Iterations", ylab="Parameter value",
          layout=c(2,4))
plotTrace(xpar$R0, axes=TRUE, div=1000)

MCMC Results for Biomass

Description

Markov chain Monte Carlo results from stock assessment of cod (Gadus morhua) in Icelandic waters, showing estimated biomass by year in tonnes.

Usage

xbio

Format

Data frame containing 1000 rows and 34 columns (years 1971 to 2004).

Details

Each column contains the results of 1 million MCMC iterations, after thinning to every 1000th iteration.

The MCMC analysis started at the best fit, so no burn-in period was discarded.

Note

Biomass is the total weight of all individuals in a population, in this case ages 4 and older.

This data frame is a subset of the xmcmc list from the scape package, which contains further documentation about the data and model. More specifically, xbio <- xmcmc$B.

The MCMC analysis was run using the AD Model Builder software (http://www.admb-project.org/).

References

Fournier, D. A., Skaug, H. J., Ancheta, J., Ianelli, J., Magnusson, A., Maunder, M. N., Nielsen, A. and Sibert, J. (2012) AD Model Builder: using automatic differentiation for statistical inference of highly parameterized complex nonlinear models. Optimization Methods and Software, 27, 233–249.

Magnusson, A., Punt, A. E. and Hilborn, R. (2013) Measuring uncertainty in fisheries stock assessment: the delta method, bootstrap, and MCMC. Fish and Fisheries, 14, 325–342.

See Also

xpar (parameters), xrec (recruitment), xbio (biomass), and xpro (projected future biomass) are MCMC data frames to explore.

plotMCMC-package gives an overview of the package.

Examples

plotDens(xbio$"2004", points=TRUE, div=1000, main="2004\n",
         xlab="Biomass age 4+ (1000 t)", tick.number=6, strip=FALSE)

plotQuant(xbio, div=1000, xlab="Year", ylab="Biomass age 4+ (kt)")
plotQuant(xbio, style="bars", div=1000, sfrac=0, xlab="Year",
          ylab="Biomass age 4+ (kt)")
plotQuant(xbio, style="lines", div=1000, xlab="Year",
          ylab="Biomass age 4+ (kt)")

MCMC Results for Model Parameters

Description

Markov chain Monte Carlo results from stock assessment of cod (Gadus morhua) in Icelandic waters, showing estimated model parameters.

Usage

xpar

Format

Data frame containing 1000 rows and 8 columns:

Details

Each column contains the results of 1 million MCMC iterations, after thinning to every 1000th iteration.

The MCMC analysis started at the best fit, so no burn-in period was discarded.

Note

This data frame is a subset of the xmcmc list from the scape package, which contains further documentation about the data and model. More specifically, xpar <- xmcmc$P.

The MCMC analysis was run using the AD Model Builder software (http://www.admb-project.org/).

References

Fournier, D. A., Skaug, H. J., Ancheta, J., Ianelli, J., Magnusson, A., Maunder, M. N., Nielsen, A. and Sibert, J. (2012) AD Model Builder: using automatic differentiation for statistical inference of highly parameterized complex nonlinear models. Optimization Methods and Software, 27, 233–249.

Magnusson, A., Punt, A. E. and Hilborn, R. (2013) Measuring uncertainty in fisheries stock assessment: the delta method, bootstrap, and MCMC. Fish and Fisheries, 14, 325–342.

See Also

xpar (parameters), xrec (recruitment), xbio (biomass), and xpro (projected future biomass) are MCMC data frames to explore.

plotMCMC-package gives an overview of the package.

Examples

plotTrace(xpar, xlab="Iterations", ylab="Parameter value",
          layout=c(2,4))
plotTrace(xpar$R0, axes=TRUE, div=1000)

plotAuto(xpar$R0)
plotAuto(xpar$R0, thin=10)
plotAuto(xpar, lag.max=50, ann=FALSE, axes=FALSE)

plotCumu(xpar$R0, main="R0")
plotCumu(xpar$cSfull, main="cSfull")
plotCumu(xpar, probs=c(0.25,0.75), ann=FALSE, axes=FALSE)

plotSplom(xpar, pch=".")

plotDens(xpar, xlab="Parameter value", ylab="Posterior density\n")

MCMC Results for Future Projections

Description

Markov chain Monte Carlo results from stock assessment of cod (Gadus morhua) in Icelandic waters, showing future projected biomass in tonnes.

Usage

xpro

Format

Data frame containing 1000 rows and 4 columns (years 2004 to 2007).

Details

Each column contains the results of 1 million MCMC iterations, after thinning to every 1000th iteration.

The MCMC analysis started at the best fit, so no burn-in period was discarded.

Note

The projections are based on a fixed harvest rate, where 25% of the biomass (ages 4 and older) is caught every year.

This data frame is a subset of the xproj list from the scape package, which contains further documentation about the data and model. More specifically, xpro <- xproj$"0.25".

The MCMC analysis was run using the AD Model Builder software (http://www.admb-project.org/).

References

Fournier, D. A., Skaug, H. J., Ancheta, J., Ianelli, J., Magnusson, A., Maunder, M. N., Nielsen, A. and Sibert, J. (2012) AD Model Builder: using automatic differentiation for statistical inference of highly parameterized complex nonlinear models. Optimization Methods and Software, 27, 233–249.

Magnusson, A., Punt, A. E. and Hilborn, R. (2013) Measuring uncertainty in fisheries stock assessment: the delta method, bootstrap, and MCMC. Fish and Fisheries, 14, 325–342.

See Also

xpar (parameters), xrec (recruitment), xbio (biomass), and xpro (projected future biomass) are MCMC data frames to explore.

plotMCMC-package gives an overview of the package.

Examples

plotQuant(xpro, axes=1:2, div=1000, xlab="Year",
          ylab="Biomass age 4+ (kt)")

plotSplom(xpro, axes=TRUE, between=1, div=1000, main="Future biomass",
          cex.labels=1.5, pch=".", cex=3)

MCMC Results for Recruitment

Description

Markov chain Monte Carlo results from stock assessment of cod (Gadus morhua) in Icelandic waters, showing estimated recruitment by year.

Usage

xrec

Format

Data frame containing 1000 rows and 33 columns (years 1970 to 2002).

Details

Each column contains the results of 1 million MCMC iterations, after thinning to every 1000th iteration.

The MCMC analysis started at the best fit, so no burn-in period was discarded.

Note

Recruitment is the size of a cohort (year class), in this case thousands of one-year-olds.

For example, xrec$"1980" is the estimated number of one-year-olds in 1981, the cohort that hatched in 1980.

This data frame is a subset of the xmcmc list from the scape package, which contains further documentation about the data and model. More specifically, xrec <- xmcmc$R.

The MCMC analysis was run using the AD Model Builder software (http://www.admb-project.org/).

References

Fournier, D. A., Skaug, H. J., Ancheta, J., Ianelli, J., Magnusson, A., Maunder, M. N., Nielsen, A. and Sibert, J. (2012) AD Model Builder: using automatic differentiation for statistical inference of highly parameterized complex nonlinear models. Optimization Methods and Software, 27, 233–249.

Magnusson, A., Punt, A. E. and Hilborn, R. (2013) Measuring uncertainty in fisheries stock assessment: the delta method, bootstrap, and MCMC. Fish and Fisheries, 14, 325–342.

See Also

xpar (parameters), xrec (recruitment), xbio (biomass), and xpro (projected future biomass) are MCMC data frames to explore.

plotMCMC-package gives an overview of the package.

Examples

plotQuant(xrec, names=substring(names(xrec),3), div=1000, xlab="Year",
          ylab="Recruitment (million one-year-olds)")