## File Name: amh_plot.R
## File Version: 0.342
#*** plot results of objects of class amh
amh_plot <- function( x, conflevel=.95, digits=3, lag.max=.1, col.smooth="red",
lwd.smooth=2, col.split="blue", lwd.split=2, lty.split=1,
col.ci="orange", cex.summ=1, ask=FALSE, ... )
{
object <- x # rename x into object
mcmcobj <- object$mcmcobj
lag.max <- round( nrow(mcmcobj) * lag.max )
round.summ <- digits
# mcmcobj <- (object$mcmcobj)[[1]]
lag.max <- min( nrow(mcmcobj), lag.max )
# index vector
a1 <- attr(mcmcobj,'mcpar')
iterindex <- seq(a1[1], a1[2], a1[3] )
# smcmcobj <- object$summary.mcmcobj
smcmcobj <- object$amh_summary
VV <- ncol(mcmcobj)
ci.quant <- - stats::qnorm( (1-conflevel)/2 )
graphics::par( mfrow=c(2,2))
for (vv in 1L:VV){
x.vv <- as.vector( mcmcobj[,vv] )
parm.vv <- colnames(mcmcobj)[vv]
sparm.vv <- smcmcobj[ smcmcobj$parameter==parm.vv, ]
#*** traceplot
graphics::plot( iterindex, x.vv, type='l',
main=paste0( 'Traceplot of ', parm.vv ),
xlab='Iterations', ylab='', ... )
x1 <- as.numeric( x.vv )
xmin <- min(x1)
xmax <- max(x1)
# 3 splits for chain
NX <- length(x.vv)
NS <- NX / 3
g1 <- round(seq( 0, NX, length=4 ))
for (ii in 1L:3){
i1 <- seq(g1[ii]+1, g1[ii+1] )
m1 <- mean( x1[ i1 ] )
graphics::lines( iterindex[ i1 ], rep( m1, length(i1) ), col=col.split,
lwd=lwd.split, lty=lty.split )
}
# include moving average here!!
l1 <- sirt_moving_average(x1, B=round( lag.max / 2 ), fill=FALSE)
graphics::lines( iterindex, l1, col=col.smooth, lwd=lwd.smooth )
#*** density estimate
graphics::plot( stats::density( x.vv ), main=paste0( 'Density of ', parm.vv ) )
c1 <- stats::quantile( x1, ( 1 - conflevel ) / 2 )
c2 <- stats::quantile( x1, 1 - ( 1 - conflevel ) / 2 )
graphics::lines( c(c1,c2), c(0,0), col=col.ci, lwd=3 )
graphics::points( sparm.vv$Mean, 0, pch=17, col=col.ci, cex=1.5)
#*** plot autocorrelation function
mtitle <- paste0( 'Autocorrelation of ', parm.vv )
m1 <- stats::acf( x.vv, lag.max=lag.max, plot=FALSE)
acf1 <- m1$acf[,1,1]
thin_lag <- round( mean( diff(iterindex) ) )
iter_vv <- seq( 0, lag.max )*thin_lag
# blue dashed line at
bd <- .05
ylim <- c( min( acf1, - bd), 1 )
graphics::plot( iter_vv, acf1, xlab='Lag', ylab='ACF',
main=mtitle, type='n', ylim=ylim)
NL <- length(iter_vv)
for (hh in 1L:NL){
graphics::lines( rep( iter_vv[hh],2), c(0, acf1[hh]) )
}
graphics::abline( h=bd, col='blue', lty=2)
graphics::abline( h=-bd, col='blue', lty=2)
#***
# numerical summary
graphics::plot( c(0,1), c(0,1), axes=FALSE, xlab='', ylab='',
main=paste0( 'Summary of ', parm.vv ), type='n', ...)
x0 <- 0
y0 <- 0
heights.summ <- c( .05, .15, .25, .35, .45, .55, .65, .75)
graphics::text( x0 + .0015, y0 + heights.summ[8], 'Posterior Mean=',
cex=cex.summ, pos=4)
graphics::text( x0 + .5, y0 + heights.summ[8],
paste0( sirt_format_numb( x=mean( x1 ), digits=round.summ) ), pos=4 )
hvv <- heights.summ[7]
graphics::text( x0 + .0015, y0 + hvv, 'Posterior Mode=', cex=cex.summ, pos=4)
graphics::text( x0 + .5, y0 + hvv,
paste0( sirt_format_numb( x=sparm.vv$MAP, digits=round.summ) ), pos=4 )
graphics::text( x0 + .0015, y0 + heights.summ[6], 'Posterior SD=',
cex=cex.summ, pos=4)
graphics::text( x0 + .5, y0 + heights.summ[6],
paste0( sirt_format_numb( x=stats::sd( x1 ), digits=round.summ) ), pos=4 )
hvv <- heights.summ[5]
graphics::text( x0 + .0015, y0 + hvv,
paste( round(100*conflevel ), '% Credibility Interval=',sep=''),
cex=cex.summ, pos=4 )
hvv <- heights.summ[4]
ci.lower <- sirt_format_numb( stats::quantile( x1, ( 1 - conflevel ) / 2 ),
digits=round.summ )
ci.upper <- sirt_format_numb( stats::quantile( x1, 1-( 1 - conflevel ) / 2 ),
digits=round.summ )
graphics::text( x0 + .25, y0 + hvv,
paste( '[', ci.lower, ',', ci.upper, ']', sep=''),
cex=cex.summ, pos=4)
hvv <- heights.summ[3]
graphics::text( x0 + .0015, y0 + hvv, 'Rhat=', cex=cex.summ, pos=4)
graphics::text( x0 + .5, y0 + hvv,
paste0( sirt_format_numb( x=sparm.vv$Rhat, digits=3) ), pos=4 )
hvv <- heights.summ[2]
graphics::text( x0 + .0015, y0 + hvv, 'SERatio=', cex=cex.summ, pos=4)
graphics::text( x0 + .5, y0 + hvv,
paste0( sirt_format_numb( x=sparm.vv$SERatio, digits=3)), pos=4 )
hvv <- heights.summ[1]
graphics::text( x0 + .0015, y0 + hvv, 'Effective Sample Size=',
cex=cex.summ, pos=4)
graphics::text( x0 + .705, y0 + hvv,
paste0( sirt_format_numb( x=sparm.vv$effSize, digits=1)), pos=4 )
graphics::par(ask=ask)
}
graphics::par(mfrow=c(1,1))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.