baujat.rma <- function(x, xlim, ylim, xlab, ylab, cex, symbol="ids", grid=TRUE, progbar=FALSE, ...) {
mstyle <- .get.mstyle()
.chkclass(class(x), must="rma", notav=c("rma.glmm", "rma.mv", "robust.rma", "rma.ls", "rma.gen", "rma.uni.selmodel"))
na.act <- getOption("na.action")
on.exit(options(na.action=na.act), add=TRUE)
if (!is.element(na.act, c("na.omit", "na.exclude", "na.fail", "na.pass")))
stop(mstyle$stop("Unknown 'na.action' specified under options()."))
if (x$k == 1L)
stop(mstyle$stop("Stopped because k = 1."))
if (is.null(x$X.f))
stop(mstyle$stop("Information needed to construct the plot is not available in the model object."))
.start.plot()
### grid argument can either be a logical or a color
if (is.logical(grid))
gridcol <- .coladj(par("bg","fg"), dark=c(0.2,-0.6), light=c(-0.2,0.6))
if (is.character(grid)) {
gridcol <- grid
grid <- TRUE
}
#########################################################################
### set up vectors to store results in
delpred <- rep(NA_real_, x$k.f)
vdelpred <- rep(NA_real_, x$k.f)
### predicted values under the full model
pred.full <- x$X.f %*% x$beta
### elements that need to be returned
outlist <- "coef.na=coef.na, beta=beta, vb=vb"
### note: skipping NA cases
### also: it is possible that model fitting fails, so that generates more NAs (these NAs will always be shown in output)
if (progbar)
pbar <- pbapply::startpb(min=0, max=x$k.f)
for (i in seq_len(x$k.f)) {
if (progbar)
pbapply::setpb(pbar, i)
if (!x$not.na[i])
next
if (inherits(x, "rma.uni"))
res <- try(suppressWarnings(.do.call(rma.uni, yi=x$yi.f, vi=x$vi.f, weights=x$weights.f, mods=x$X.f, intercept=FALSE, method=x$method, weighted=x$weighted, test=x$test, level=x$level, tau2=ifelse(x$tau2.fix, x$tau2, NA), control=x$control, subset=-i, skipr2=TRUE, outlist=outlist)), silent=TRUE)
if (inherits(x, "rma.mh")) {
if (is.element(x$measure, c("RR","OR","RD"))) {
res <- try(suppressWarnings(.do.call(rma.mh, ai=x$outdat.f$ai, bi=x$outdat.f$bi, ci=x$outdat.f$ci, di=x$outdat.f$di, measure=x$measure, add=x$add, to=x$to, drop00=x$drop00, correct=x$correct, level=x$level, subset=-i, outlist=outlist)), silent=TRUE)
} else {
res <- try(suppressWarnings(.do.call(rma.mh, x1i=x$outdat.f$x1i, x2i=x$outdat.f$x2i, t1i=x$outdat.f$t1i, t2i=x$outdat.f$t2i, measure=x$measure, add=x$add, to=x$to, drop00=x$drop00, correct=x$correct, level=x$level, subset=-i, outlist=outlist)), silent=TRUE)
}
}
if (inherits(x, "rma.peto"))
res <- try(suppressWarnings(.do.call(rma.peto, ai=x$outdat.f$ai, bi=x$outdat.f$bi, ci=x$outdat.f$ci, di=x$outdat.f$di, add=x$add, to=x$to, drop00=x$drop00, level=x$level, subset=-i)), silent=TRUE)
if (inherits(res, "try-error"))
next
### removing an observation could lead to a model coefficient becoming inestimable (for 'rma.uni' objects)
if (any(res$coef.na))
next
Xi <- matrix(x$X.f[i,], nrow=1)
delpred[i] <- Xi %*% res$beta
vdelpred[i] <- Xi %*% tcrossprod(res$vb,Xi)
}
if (progbar)
pbapply::closepb(pbar)
yhati <- (delpred - pred.full)^2 / vdelpred
#########################################################################
### x-axis values (use 'na.pass' to make sure we get a vector of length k.f)
options(na.action = "na.pass")
xhati <- resid(x)^2 / (x$tau2.f + x$vi.f)
options(na.action = na.act)
#########################################################################
### set some defaults (if not specified)
if (missing(cex))
cex <- par("cex") * 0.8
if (missing(xlab)) {
if (is.element(x$method, c("FE","EE","CE"))) {
xlab <- ifelse(x$int.only, "Contribution to Overall Heterogeneity", "Contribution to Residual Heterogeneity")
} else {
xlab <- "Squared Pearson Residual"
}
}
if (missing(ylab))
ylab <- ifelse(x$int.only, "Influence on Overall Result", "Influence on Fitted Value")
if (missing(xlim))
xlim <- range(xhati, na.rm=TRUE)
if (missing(ylim))
ylim <- range(yhati, na.rm=TRUE)
#########################################################################
### draw empty plot
plot(NA, xlab=xlab, ylab=ylab, xlim=xlim, ylim=ylim, ...)
### add grid (and redraw box)
if (.isTRUE(grid)) {
grid(col=gridcol)
box(...)
}
if (is.numeric(symbol)) {
symbol <- .expand1(symbol, x$k.all)
if (length(symbol) != x$k.all)
stop(mstyle$stop(paste0("Length of the 'symbol' argument (", length(symbol), ") does not correspond to the size of the original dataset (", x$k.all, ").")))
symbol <- .getsubset(symbol, x$subset)
points(x=xhati, y=yhati, cex=cex, pch=symbol, ...)
}
if (is.character(symbol) && symbol=="ids")
text(xhati, yhati, x$ids, cex=cex, ...)
if (is.character(symbol) && symbol=="slab")
text(xhati, yhati, x$slab, cex=cex, ...)
#########################################################################
sav <- data.frame(x=xhati[x$not.na], y=yhati[x$not.na], ids=x$ids[x$not.na], slab=x$slab[x$not.na], stringsAsFactors=FALSE)
invisible(sav)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.