Nothing
plot.mitml <- function(x, print = c("beta", "beta2", "psi", "sigma"), pos = NULL, group = "all",
trace = c("imputation", "burnin", "all"), thin = 1, smooth = 3,
n.Rhat = 3, export = c("none", "png", "pdf"), dev.args = list(),
...){
# plot method for objects of class "mitml"
# retrieve data and variable names (predictors)
vrs <- x$model
clus <- x$model$clus
pvrs <- seq_along(attr(vrs, "full.names")$pvrs)
qvrs <- seq_along(attr(vrs, "full.names")$qvrs)
names(pvrs) <- attr(vrs, "full.names")$pvrs
names(qvrs) <- attr(vrs, "full.names")$qvrs
isML <- attr(x$model, "is.ML")
isL2 <- attr(x$model, "is.L2")
if(isL2){
pvrs.L2 <- seq_along(attr(vrs, "full.names")$pvrs.L2)
names(pvrs.L2) <- attr(vrs, "full.names")$pvrs.L2
}
# match arguments
print <- match.arg(print, several.ok = TRUE)
trace <- match.arg(trace)
export <- match.arg(export)
# check for random L1
rl1 <- x$random.L1 == "full"
# parameter chains (for backwards compatibility)
kc <- x$keep.chains
if(is.null(kc)) kc <- "full"
# check print and position for selected parameters
if(!is.null(pos) & length(print)>1){
pos <- NULL
warning("The 'pos' argument may only be used when 'print' is cleary defined as one of 'beta', 'beta2', 'psi', or 'sigma' (see '?plot').")
}
# grouping
grp.labels <- unique(attr(x$data, "group"))
if(is.numeric(group)) grp.labels <- grp.labels[group]
grp <- length(grp.labels)
# export, graphical parameters
if(export != "none"){
wd <- getwd()
out <- file.path(wd, "mitmlPlots")
if(!file.exists(out)) dir.create(out)
}else{
do.call(dev.new, dev.args)
devAskNewPage(ask = FALSE)
}
oldpar <- par(no.readonly = TRUE)
# ***
# start plotting
#
for(gg in 1:grp){
# grouping
if(grp>1){
glab <- paste(",Group:", grp.labels[gg], sep = "")
gfile <- paste("Group-", grp.labels[gg], "_", sep = "")
}else{
glab <- gfile <- ""
}
# expand targets for multiple categories
yvrs <- vrs$yvrs
yvrs.L2 <- vrs$yvrs.L2
# ... level 1
cvrs <- intersect(yvrs, attr(x$data, "cvrs"))
nc <- length(cvrs)
if(length(cvrs)>=1){
yvrs <- c(yvrs[!yvrs %in% cvrs], cvrs)
for(cc in 1:nc){
cv <- cvrs[cc]
ci <- which(yvrs == cv)
yi <- 1:length(yvrs)
nlev <- attr(x$data, "levels")[gg, cc]
if(nlev>2){
newy <- paste0(cv, 1:(nlev-1))
}else{
newy <- cv
}
sel0 <- yi[yi<ci]
sel1 <- yi[yi>ci]
yvrs <- c(yvrs[sel0], newy, yvrs[sel1])
}
}
ynam <- yvrs
yvrs <- seq_along(yvrs)
names(yvrs) <- ynam
# ... level 2
if(isL2){
cvrs.L2 <- intersect(yvrs.L2, attr(x$data, "cvrs"))
nc.L2 <- length(cvrs.L2)
if(length(cvrs.L2)>=1){
yvrs.L2 <- c(yvrs.L2[!yvrs.L2 %in% cvrs.L2], cvrs.L2)
for(cc in 1:nc.L2){
cv <- cvrs.L2[cc]
ci <- which(yvrs.L2 == cv)
yi <- 1:length(yvrs.L2)
nlev <- attr(x$data, "levels")[gg, nc+cc]
if(nlev>2){
newy <- paste0(cv, 1:(nlev-1))
}else{
newy <- cv
}
sel0 <- yi[yi<ci]
sel1 <- yi[yi>ci]
yvrs.L2 <- c(yvrs.L2[sel0], newy, yvrs.L2[sel1])
}
}
ynam <- yvrs.L2
yvrs.L2 <- seq_along(yvrs.L2)
names(yvrs.L2) <- ynam
}
# number of iterations
n <- dim(x$par.burnin[["beta"]])[3]+dim(x$par.imputation[["beta"]])[3]
nb <- dim(x$par.burnin[["beta"]])[3]
ni <- dim(x$par.imputation[["beta"]])[3]
niter <- x$iter[["iter"]]
# thinned-sample indicators
s <- seq.int(thin, n, by = thin)
sb <- seq.int(thin, nb, by = thin)
si <- seq.int(thin, ni, by = thin)
lag <- ceiling(niter/thin)
# *** plots for fixed regression coefficients at level 1
#
if("beta" %in% print){
# check if pos is badly defined
if(!is.null(pos)){
if(pos[1] > max(pvrs) | pos[1] < min(pvrs) | pos[2] > max(yvrs) | pos[2] < min(yvrs)){
.restoreDevice(oldpar, export, close = TRUE)
stop("There is no entry [", pos[1], ",", pos[2], "] in 'beta'.")
}
}
for(ic in yvrs){
for(ir in pvrs){
# skip if individual parameters requested
if(!is.null(pos)){
if(!(pos[1] == ir & pos[2] == ic)) next
}
if(export != "none"){
filename <- paste("BETA_", gfile, names(yvrs[ic]), "_ON_", names(pvrs[ir]), ".", export, sep = "")
filename <- gsub("[(),]", "", filename)
filename <- gsub("[[:space:]]", "-", filename)
out.args <- c(list(file = file.path(out, filename)), dev.args)
do.call(export, out.args)
}
layout(matrix(c(1, 2, 3, 4), 2, 2), c(5, 1), c(1.13, 1))
# choose section of trace
switch(trace,
imputation={
trc <- x$par.imputation[["beta"]][ir,ic,,gg][si]
},
burnin={
trc <- x$par.burnin[["beta"]][ir,ic,,gg][sb]
},
all={
trc <- c(x$par.burnin[["beta"]][ir,ic,,gg],
x$par.imputation[["beta"]][ir,ic,,gg])[s]
}
)
# trace plot
par(mar = c(3, 3, 2, 0)+0.5, mgp = c(2, 1, 0), font.lab = 2)
ymin <- min(trc)
ymax <- max(trc)
yr <- ymax-ymin
plot(trc, type = ifelse(trace == "all", "n", "l"), ylab = "Trace", xlab = "Iteration",
xaxt = "n", ylim = c(ymin-yr*.03, ymax+yr*.03),
...)
if(trace == "all"){
lines(which(s<=nb), trc[s<=nb], col = "grey75")
lines(which(s>=nb), trc[s>=nb], col = "black")
}
axt <- axTicks(1)
title(main = paste("Beta [", ir, ",", ic, glab, "]: ", names(yvrs[ic]), " ON ", names(pvrs[ir]), sep = ""), cex.main = 1)
if(trace == "imputation"){
axl <- sprintf("%d", thin*(axt+nb))
}else{
axl <- sprintf("%d", thin*axt)
}
axis(side = 1, at = axt, labels = axl)
# trend line for trace (moving window average)
if(all(is.numeric(smooth), smooth>0)){
B <- floor(niter/(smooth*thin))
mwa <- .movingAverage(trc, B, fill = TRUE)
lines(mwa, col = "grey60")
}
# blue line
if(trace == "all") abline(v = ceiling(nb/thin), col = "blue")
# further plots
if(trace == "burnin"){
drw <- x$par.burnin[["beta"]][ir, ic, sb, gg]
}else{
drw <- x$par.imputation[["beta"]][ir, ic, si, gg]
}
# autocorrelation plot
par(mar = c(3, 3, 1, 0)+0.5)
ac <- acf(drw, lag.max = lag+2, plot = F)
plot(ac[1:lag], ylim = c(-.1, 1), yaxt = "n", main = NULL, ylab = "ACF", ci = 0,
...)
axis(side = 2, at = c(0, .5, 1))
abline(h = c(-.1, .1), col = "blue")
# kernel density plot
par(mar = c(3, 0, 2, 0)+0.5, mgp = c(2, 0, 0))
ddrw <- density(drw)
plot(x = ddrw$y, y = ddrw$x, type = "l", xaxt = "n", yaxt = "n", xlab = "", ylab = "",
ylim = c(ymin-yr*.03, ymax+yr*.03),
...)
# posterior summary
par(mar = c(1, -0.5, 0, -0.5)+0.5)
plot.new()
text(0, 0.5, paste("EAP: ", sprintf(fmt = "%.3f", mean(drw)), "\n",
"MAP: ", sprintf(fmt = "%.3f", ddrw$x[which.max(ddrw$y)]), "\n",
"SD: ", sprintf(fmt = "%.3f", sd(drw)), "\n",
"2.5%: ", sprintf(fmt = "%.3f", quantile(drw, .025)), "\n",
"97.5%: ", sprintf(fmt = "%.3f", quantile(drw, .975)), "\n",
"Rhat: ", sprintf(fmt = "%.3f", .GelmanRubin(t(drw), n.Rhat)), "\n",
"ACF-k: ", sprintf(fmt = "%.3f", .smoothedACF(ac, k = lag, sd=.5)), "\n",
sep = ""), adj = c(0, .5), cex=.8, family = "mono", font = 2, ...)
if(export != "none"){
dev.off()
}else{
devAskNewPage(ask = TRUE)
}
}}}
# *** plots for fixed regression coefficients at level 2
#
if(isL2 & "beta2" %in% print){
# check if pos is badly defined
if(!is.null(pos)){
if(pos[1] > max(pvrs.L2) | pos[1] < min(pvrs.L2) | pos[2] > max(yvrs.L2) | pos[2] < min(yvrs.L2)){
.restoreDevice(oldpar, export, close = TRUE)
stop("There is no entry [", pos[1], ",", pos[2], "] in 'beta2'.")
}
}
for(ic in yvrs.L2){
for(ir in pvrs.L2){
# skip if individual parameters requested
if(!is.null(pos)){
if(!(pos[1] == ir & pos[2] == ic)) next
}
if(export != "none"){
filename <- paste("BETA2_", gfile, names(yvrs.L2[ic]), "_ON_", names(pvrs.L2[ir]), ".", export, sep = "")
filename <- gsub("[(),]", "", filename)
filename <- gsub("[[:space:]]", "-", filename)
out.args <- c(list(file = file.path(out, filename)), dev.args)
do.call(export, out.args)
}
layout(matrix(c(1, 2, 3, 4), 2, 2), c(5, 1), c(1.13, 1))
# choose section of trace
switch(trace,
imputation={
trc <- x$par.imputation[["beta2"]][ir,ic,,gg][si]
},
burnin={
trc <- x$par.burnin[["beta2"]][ir,ic,,gg][sb]
},
all={
trc <- c(x$par.burnin[["beta2"]][ir,ic,,gg],
x$par.imputation[["beta2"]][ir,ic,,gg])[s]
}
)
# trace plot
par(mar = c(3, 3, 2, 0)+0.5, mgp = c(2, 1, 0), font.lab = 2)
ymin <- min(trc)
ymax <- max(trc)
yr <- ymax-ymin
plot(trc, type = ifelse(trace == "all", "n", "l"), ylab = "Trace", xlab = "Iteration",
xaxt = "n", ylim = c(ymin-yr*.03, ymax+yr*.03),
...)
if(trace == "all"){
lines(which(s<=nb), trc[s<=nb], col = "grey75")
lines(which(s>=nb), trc[s>=nb], col = "black")
}
axt <- axTicks(1)
title(main = paste("Beta2 [", ir, ",", ic, glab, "]: ", names(yvrs.L2[ic]), " ON ", names(pvrs.L2[ir]), sep = ""), cex.main = 1)
if(trace == "imputation"){
axl <- sprintf("%d", thin*(axt+nb))
}else{
axl <- sprintf("%d", thin*axt)
}
axis(side = 1, at = axt, labels = axl)
# trend line for trace (moving window average)
if(all(is.numeric(smooth), smooth>0)){
B <- floor(niter/(smooth*thin))
mwa <- .movingAverage(trc, B, fill = TRUE)
lines(mwa, col = "grey60")
}
# blue line
if(trace == "all") abline(v = ceiling(nb/thin), col = "blue")
# further plots
if(trace == "burnin"){
drw <- x$par.burnin[["beta2"]][ir,ic,sb,gg]
}else{
drw <- x$par.imputation[["beta2"]][ir,ic,si,gg]
}
# autocorrelation plot
par(mar = c(3, 3, 1, 0)+0.5)
ac <- acf(drw, lag.max = lag+2, plot = F)
plot(ac[1:lag], ylim = c(-.1, 1), yaxt = "n", main = NULL, ylab = "ACF", ci = 0,
...)
axis(side = 2, at = c(0, .5, 1))
abline(h = c(-.1, .1), col = "blue")
# kernel density plot
par(mar = c(3, 0, 2, 0)+0.5, mgp = c(2, 0, 0))
ddrw <- density(drw)
plot(x = ddrw$y, y = ddrw$x, type = "l", xaxt = "n", yaxt = "n", xlab = "", ylab = "",
ylim = c(ymin-yr*.03, ymax+yr*.03),
...)
# posterior summary
par(mar = c(1, -0.5, 0, -0.5)+0.5)
plot.new()
text(0, 0.5, paste("EAP: ", sprintf(fmt = "%.3f", mean(drw)), "\n",
"MAP: ", sprintf(fmt = "%.3f", ddrw$x[which.max(ddrw$y)]), "\n",
"SD: ", sprintf(fmt = "%.3f", sd(drw)), "\n",
"2.5%: ", sprintf(fmt = "%.3f", quantile(drw, .025)), "\n",
"97.5%: ", sprintf(fmt = "%.3f", quantile(drw, .975)), "\n",
"Rhat: ", sprintf(fmt = "%.3f", .GelmanRubin(t(drw), n.Rhat)), "\n",
"ACF-k: ", sprintf(fmt = "%.3f", .smoothedACF(ac, k = lag, sd=.5)), "\n",
sep = ""), adj = c(0, .5), cex=.8, family = "mono", font = 2, ...)
if(export != "none"){
dev.off()
}else{
devAskNewPage(ask = TRUE)
}
}}}
# *** plots for random effects' variance components
#
if(isML & "psi" %in% print){
# joint set of variables at level 1 and 2
yvrs.comb <- c(yvrs, if(isL2) yvrs.L2+length(yvrs))
# index matrix
bvec <- t(expand.grid(qvrs, yvrs))
if(isL2) bvec <- cbind(bvec, t(expand.grid(1, yvrs.L2+length(yvrs))))
# attempt to fix pos if badly defined
if(!is.null(pos)){
pos0 <- pos
if(pos[2]>pos[1]){ # fix if pos is redundant/transposed
pos[1] <- pos0[2]
pos[2] <- pos0[1]
}
if(any(pos0 > max(yvrs.comb)) | any(pos0 < min(yvrs.comb))){
.restoreDevice(oldpar, export, close = TRUE)
stop("There is no entry [", pos0[1], ",", pos0[2], "] in 'psi'.")
}
if(!identical(pos, pos0))
warning("Could not use entry [", pos0[1], ",", pos0[2], "] in 'psi'. Used [", pos[1], ",", pos[2], "] instead.")
}
dpsi <- length(yvrs)*length(qvrs)
if(isL2) dpsi <- dpsi+length(yvrs.L2)
# if only "diagonal" entries, fix max. column index to 1
cpsi <- if(kc == "diagonal") 1 else dpsi
for(ic in 1:cpsi){
for(ir in ic:dpsi){
# skip if different individual parameters requested
if(!is.null(pos)){
if(!(pos[1] == ir & pos[2] == ic)) next
}
# if only "diagonal" entries, use ir for all labels
ic2 <- if(kc == "diagonal") ir else ic
# check for residual at L2
icL2 <- ic > (length(yvrs)*length(qvrs))
irL2 <- ir > (length(yvrs)*length(qvrs))
if(export != "none"){
filename <- paste0("PSI_", gfile,
names(yvrs.comb[bvec[2, ir]]),
if(!irL2) paste0("_ON_", names(qvrs[bvec[1, ir]])),
"_WITH_",
names(yvrs.comb[bvec[2, ic2]]),
if(!icL2) paste0("_ON_", names(qvrs[bvec[1, ic2]])),
".", export)
filename <- gsub("[(),]", "", filename)
filename <- gsub("[[:space:]]", "-", filename)
out.args <- c(list(file = file.path(out, filename)), dev.args)
do.call(export, out.args)
}
layout(matrix(c(1, 2, 3, 4), 2, 2), c(5, 1), c(1.13, 1))
switch(trace,
imputation={
trc <- x$par.imputation[["psi"]][ir,ic,,gg][si]
},
burnin={
trc <- x$par.burnin[["psi"]][ir,ic,,gg][sb]
},
all={
trc <- c(x$par.burnin[["psi"]][ir,ic,,gg],
x$par.imputation[["psi"]][ir,ic,,gg])[s]
}
)
# trace plot
par(mar = c(3, 3, 2, 0)+0.5, mgp = c(2, 1, 0), font.lab = 2)
ymin <- min(trc)
ymax <- max(trc)
yr <- ymax-ymin
plot(trc, type = ifelse(trace == "all", "n", "l"), ylab = "Trace", xlab = "Iteration",
xaxt = "n", ylim = c(ymin-yr*.03, ymax+yr*.03),
...)
if(trace == "all"){
lines(which(s<=nb), trc[s<=nb], col = "grey75")
lines(which(s>=nb), trc[s>=nb], col = "black")
}
title(main = paste0("Psi [", ir, ",", ic, glab, "]: ",
if(!irL2) "(",
names(yvrs.comb[bvec[2, ir]]),
if(!irL2) paste0(" ON ", names(qvrs[bvec[1, ir]]), ")"),
" WITH ",
if(!icL2) "(",
names(yvrs.comb[bvec[2, ic2]]),
if(!icL2) paste0(" ON ", names(qvrs[bvec[1, ic2]]), ")")
),
cex.main = 1)
axt <- axTicks(1)
if(trace == "imputation"){
axl <- sprintf("%d", thin*(axt+nb))
}else{
axl <- sprintf("%d", thin*axt)
}
axis(side = 1, at = axt, labels = axl)
# trend line for trace (moving window average)
if(all(is.numeric(smooth), smooth>0)){
B <- floor(niter/(smooth*thin))
mwa <- .movingAverage(trc, B, fill = TRUE)
lines(mwa, col = "grey60")
}
# blue line
if(trace == "all") abline(v = ceiling(nb/thin), col = "blue")
# further plots
if(trace == "burnin"){
drw <- x$par.burnin[["psi"]][ir,ic,sb,gg]
}else{
drw <- x$par.imputation[["psi"]][ir,ic,si,gg]
}
# autocorrelation plot
par(mar = c(3, 3, 1, 0)+0.5)
ac <- acf(drw, lag.max = lag+2, plot = F)
plot(ac[1:lag], ylim = c(-.1, 1), yaxt = "n", main = NULL, ylab = "ACF", ci = 0,
...)
axis(side = 2, at = c(0, .5, 1))
abline(h = c(-.1, .1), col = "blue")
# kernel density plot
par(mar = c(3, 0, 2, 0)+0.5, mgp = c(2, 0, 0))
ddrw <- density(drw)
plot(x = ddrw$y, y = ddrw$x, type = "l", xaxt = "n", yaxt = "n", xlab = "", ylab = "", ylim = c(ymin-yr*.03, ymax+yr*.03), ...)
# posterior summary
par(mar = c(1, -0.5, 0, -0.5)+0.5)
plot.new()
text(0, 0.5, paste("EAP: ", sprintf(fmt = "%.3f", mean(drw)), "\n",
"MAP: ", sprintf(fmt = "%.3f", ddrw$x[which.max(ddrw$y)]), "\n",
"SD: ", sprintf(fmt = "%.3f", sd(drw)), "\n",
"2.5%: ", sprintf(fmt = "%.3f", quantile(drw, .025)), "\n",
"97.5%: ", sprintf(fmt = "%.3f", quantile(drw, .975)), "\n",
"Rhat: ", sprintf(fmt = "%.3f", .GelmanRubin(t(drw), n.Rhat)), "\n",
"ACF-k: ", sprintf(fmt = "%.3f", .smoothedACF(ac, k = lag, sd=.5)), "\n",
sep = ""), adj = c(0, .5), cex=.8, family = "mono", font = 2, ...)
if(export != "none"){
dev.off()
}else{
devAskNewPage(ask = TRUE)
}
}}}
# *** plots for residual variance components
#
if("sigma" %in% print){
# cluster-specific covariance matrices stacked in rows
gind <- attr(x$data, "group") == grp.labels[gg]
clus2 <- unique(x$data[gind, clus])
clus3 <- if(rl1) seq_along(clus2) else 1
# attempt to fix pos if badly defined
if(!is.null(pos)){
pos0 <- pos
dims <- dim(x$par.imputation$sigma)
if(pos[2] > length(yvrs)){ # fix if pos is transposed
pos[1] <- pos0[2]
pos[2] <- pos0[1]
pos0 <- pos
}
if(pos[2] > ((pos[1]-1)%%length(yvrs))+1){ # fix if pos is redundant
pos[1] <- pos0[1] - pos0[1]%%length(yvrs) + pos0[2]
pos[2] <- pos0[1]%%length(yvrs)
}
if(all(pos0 > max(yvrs)) | any(pos0 < min(yvrs)) | max(pos0) > dims[1]){
.restoreDevice(oldpar, export, close = TRUE)
stop("There is no entry [", pos0[1], ",", pos0[2], "] in 'sigma'.")
}
if(!identical(pos, pos0))
warning("Could not use entry [", pos0[1], ",", pos0[2], "] in 'sigma'. Used [", pos[1], ",", pos[2], "] instead.")
}
# if only "diagonal" entries, fix max. column index to 1
csig <- if(kc == "diagonal") 1 else length(yvrs)
for(icl in clus3){
for(ic in 1:csig){
for(ir in ic:length(yvrs)){
# adjust row index for cluster-specific covariance matrices
ir2 <- ir+(icl-1)*length(yvrs)
# if only "diagonal" entries, use ir for all labels
ic2 <- if(kc == "diagonal") ir else ic
# skip if individual parameters requested
if(!is.null(pos)){
if(!(pos[1] == ir2 & pos[2] == ic)) next
}
if(export != "none"){
filename <- paste0("SIGMA_", gfile,
names(yvrs[ir]),
"_WITH_",
names(yvrs[ic2]),
if(rl1) paste0("_", clus, clus2[icl]),
".", export)
filename <- gsub("[(),]", "", filename)
filename <- gsub("[[:space:]]", "-", filename)
out.args <- c(list(file = file.path(out, filename)), dev.args)
do.call(export, out.args)
}
layout(matrix(c(1, 2, 3, 4), 2, 2), c(5, 1), c(1.13, 1))
switch(trace,
imputation={
trc <- x$par.imputation[["sigma"]][ir2,ic,,gg][si]
},
burnin={
trc <- x$par.burnin[["sigma"]][ir2,ic,,gg][sb]
},
all={
trc <- c(x$par.burnin[["sigma"]][ir2,ic,,gg],
x$par.imputation[["sigma"]][ir2,ic,,gg])[s]
}
)
# trace plots
par(mar = c(3, 3, 2, 0)+0.5, mgp = c(2, 1, 0), font.lab = 2)
ymin <- min(trc)
ymax <- max(trc)
yr <- ymax-ymin
plot(trc, type = ifelse(trace == "all", "n", "l"), ylab = "Trace", xlab = "Iteration",
xaxt = "n", ylim = c(ymin-yr*.03, ymax+yr*.03),
...)
if(trace == "all"){
lines(which(s<=nb), trc[s<=nb], col = "grey75")
lines(which(s>=nb), trc[s>=nb], col = "black")
}
title(main = paste0("Sigma [", ir2, ",", ic, glab, "]: ",
names(yvrs[ir]),
" WITH ",
names(yvrs[ic2]),
if(rl1) paste0(" [", clus, ":", clus2[icl], "]")),
cex.main = 1)
axt <- axTicks(1)
if(trace == "imputation"){
axl <- sprintf("%d", thin*(axt+nb))
}else{
axl <- sprintf("%d", thin*axt)
}
axis(side = 1, at = axt, labels = axl)
# trend line for trace (moving window average)
if(all(is.numeric(smooth), smooth>0)){
B <- floor(niter/(smooth*thin))
mwa <- .movingAverage(trc, B, fill = TRUE)
lines(mwa, col = "grey60")
}
# blue line
if(trace == "all") abline(v = ceiling(nb/thin), col = "blue")
# further plots
if(trace == "burnin"){
drw <- x$par.burnin[["sigma"]][ir2,ic,sb,gg]
}else{
drw <- x$par.imputation[["sigma"]][ir2,ic,si,gg]
}
# autocorrelation plot
par(mar = c(3, 3, 1, 0)+0.5)
ac <- acf(drw, lag.max = lag+2, plot = F)
plot(ac[1:lag], ylim = c(-.1, 1), yaxt = "n", main = NULL, ylab = "ACF", ci = 0,
...)
axis(side = 2, at = c(0, .5, 1))
abline(h = c(-.1, .1), col = "blue")
# kernel density plot
par(mar = c(3, 0, 2, 0)+0.5, mgp = c(2, 0, 0))
ddrw <- density(drw)
plot(x = ddrw$y, y = ddrw$x, type = "l", xaxt = "n", yaxt = "n", xlab = "", ylab = "", ylim = c(ymin-yr*.03, ymax+yr*.03), ...)
# posterior summary
par(mar = c(1, -0.5, 0, -0.5)+0.5)
plot.new()
text(0, 0.5, paste("EAP: ", sprintf(fmt = "%.3f", mean(drw)), "\n",
"MAP: ", sprintf(fmt = "%.3f", ddrw$x[which.max(ddrw$y)]), "\n",
"SD: ", sprintf(fmt = "%.3f", sd(drw)), "\n",
"2.5%: ", sprintf(fmt = "%.3f", quantile(drw, .025)), "\n",
"97.5%: ", sprintf(fmt = "%.3f", quantile(drw, .975)), "\n",
"Rhat: ", sprintf(fmt = "%.3f", .GelmanRubin(t(drw), n.Rhat)), "\n",
"ACF-k: ", sprintf(fmt = "%.3f", .smoothedACF(ac, k = lag, sd=.5)), "\n",
sep = ""), adj = c(0, .5), cex=.8, family = "mono", font = 2, ...)
if(export != "none"){
dev.off()
}else{
devAskNewPage(ask = TRUE)
}
}}}}
}
plot.new()
par(oldpar)
if(export == "none") devAskNewPage(ask = FALSE)
dev.off()
invisible()
}
# restore and shut down parameters upon error
.restoreDevice <- function(pars, export, close = TRUE){
par(pars)
if(export == "none") devAskNewPage(ask = FALSE)
if(close) dev.off()
invisible()
}
# moving window average for time series
.movingAverage <- function(x, B, fill = TRUE){
x1 <- cumsum(x)
N <- length(x)
y <- rep(NA, N)
i <- seq(B+1 , N-B)
xdiff <- x1[ -seq(1, B) ] - x1[ -seq(N-B+1, N) ]
xdiff <- xdiff[ - seq(1, B) ]
y[i] <- ( x1[i] + xdiff - c(0, x1[ -seq(N-2*B, N) ]) ) / (2*B+1)
# fill NAs at beginning and end of time series
if(fill){
j <- seq(0, B-1)
ybeg <- sapply(j, function(z) sum( x[ seq(1, (2*z+1)) ]) / (2*z+1) )
yend <- sapply(rev(j), function(z) sum( x[ seq(N-2*z, N) ] ) / (2*z+1) )
y[j+1] <- ybeg
y[rev(N-j)] <- yend
}
y
}
# lag-k autocorrelation smoothed by values of a normal density
.smoothedACF <- function(x, k, sd=.5){
x0 <- x$ac[-1, 1, 1]
n <- length(x0)
add <- n-k
x0 <- x0[(k-add):n]
# weights based on normal density
w <- dnorm(-add:add, 0, sd)
y <- sum( x0 * (w/sum(w)) )
y
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.