Predict method for Robust Bayesian Meta-Analysis Fits

predict.RoBMA predicts values based on the RoBMA model. Only available for normal-normal models estimated using the spike-and-slab algorithm (i.e., algorithm = "ss").

Usage

# S3 method for class 'RoBMA'
predict(
  object,
  newdata = NULL,
  type = "response",
  conditional = FALSE,
  output_scale = NULL,
  probs = c(0.025, 0.975),
  incorporate_publication_bias = TRUE,
  as_samples = FALSE,
  ...
)

Arguments

object: a fitted RoBMA object
newdata: a data.frame (if prediction for a meta-regression is performed) or a list named list with the effect size measure and variability metrics (if prediction for a meta-analysis is performed) for new studies. Note that the input has to corresponds to the format and naming that was used to estimate the original fit. Defaults to NULL which corresponds to prediction for the observed data.
type: type of prediction to be performed. Defaults to "response" which produces predictions for the observed effect size estimates. Alternatives are "terms" which produces the mean effect size estimate at the given predictors levels (not accounting for the random-effects) and "effect" which predicts the distribution of the true study effects at the given predictors levels (i.e., incorporating heterogeneity into "terms").
conditional: show the conditional estimates (assuming that the alternative is true). Defaults to FALSE. Only available for type == "ensemble".
output_scale: transform the meta-analytic estimates to a different scale. Defaults to NULL which returns the same scale as the model was estimated on.
probs: quantiles of the posterior samples to be displayed. Defaults to c(.025, .975)
incorporate_publication_bias: whether sampling of new values should incorporate the estimated publication bias (note that selection models do not affect the mean paramater when "terms" (equal mean parameter under normal vs. weighted likelihood equals different expectation).
as_samples: whether posterior samples instead of a summary table should be returned. Defaults to FALSE.
...: additional arguments

Value

pooled_effect returns a list of tables of class 'BayesTools_table'.

Details

Note that in contrast to predict, the type = "response" produces predictions for the new effect size estimates (instead of the true study effects). To obtain results corresponding to the metafor's predict function, use the type = "terms" to obtain the mean effect size estimate in its credible interval and type = "effect" to obtain the distribution of the true study effects (i.e., prediction interval).

The conditional estimate is calculated conditional on the presence of the effect (in meta-analysis) or the intercept (in meta-regression).

Examples

if (FALSE) { # \dontrun{
require(metafor)
dat <- escalc(measure = "OR", ai = tpos, bi = tneg, ci = cpos, di = cneg,
              data = dat.bcg)

# fit meta-regression
robma_dat <- data.frame(
  logOR  = dat$yi,
  se     = sqrt(dat$vi),
  ablat  = dat$ablat,
  alloc  = dat$alloc
)

fit <- NoBMA.reg(~ ablat + alloc, data = robma_dat,
                 seed = 1, algorithm = "ss", parallel = TRUE)

# prediction for the mean effect, prediction interval, and posterior predictive
mean_effect          <- predict(fit, type = "terms")
prediction_interval  <- predict(fit, type = "effect")
posterior_predictive <- predict(fit, type = "response")

# visualize the estimates vs predictions
plot(NA, type = "n", xlim = c(-3, 3), ylim = c(0, nrow(robma_dat) + 1),
     xlab = "logOR", ylab = "Observation", las = 1)
points(robma_dat$logOR, seq_along(robma_dat$logOR), cex = 2, pch = 16)
points(mean_effect$estimates$Mean, seq_along(robma_dat$logOR) + 0.2,
       cex = 1.5, pch = 16, col = "blue")
sapply(seq_along(robma_dat$logOR), function(i){
  lines(c(robma_dat$logOR[i] - 1.96 * robma_dat$se[i],
          robma_dat$logOR[i] + 1.96 * robma_dat$se[i]),
        c(i, i), lwd = 2)
  lines(c(mean_effect$estimates[i,"0.025"],
          mean_effect$estimates[i,"0.975"]),
        c(i + 0.2, i + 0.2), lwd = 2, col = "blue")
  lines(c(prediction_interval$estimates[i,"0.025"],
          prediction_interval$estimates[i,"0.975"]),
        c(i + 0.3, i + 0.3), lwd = 2, lty = 2, col = "blue")
  lines(c(posterior_predictive$estimates[i,"0.025"],
          posterior_predictive$estimates[i,"0.975"]),
        c(i + 0.4, i + 0.4), lwd = 2, lty = 3, col = "blue")
})
legend("bottomright", col = c("black", rep("blue", 3)),
       lwd = 2, lty = c(1, 1, 2, 3), bty = "n",
       legend = c("Observed + CI", "Predicted + CI",
                  "Prediction Int.", "Sampling Int.")
      )

 # prediction across a lattitude
 fit2 <- NoBMA.reg(~ ablat, data = robma_dat,
                   seed = 1, algorithm = "ss", parallel = TRUE)

 new_df <- data.frame(
   logOR  = 0,   # only relevant for "response" (not plotted here)
   se     = 0.1, # only relevant for "response" (not plotted here)
   ablat  = 10:60
 )
 new_mean_effect          <- predict(fit2, newdata = new_df, type = "terms")
 new_prediction_interval  <- predict(fit2, newdata = new_df, type = "effect")

 # create bubble plot
 plot(robma_dat$ablat, robma_dat$logOR,  ylim = c(-2, 1),
      xlab = "Latitude", ylab = "logOR", las = 1, cex = 0.3/robma_dat$se, pch = 16)
 polygon(c(10:60, rev(10:60)),
         c(new_mean_effect$estimates[,"0.025"],
           rev(new_mean_effect$estimates[,"0.975"])),
         col = rgb(0, 0, 1, alpha = 0.2), border = NA)
 polygon(c(10:60, rev(10:60)),
         c(new_prediction_interval$estimates[,"0.025"],
           rev(new_prediction_interval$estimates[,"0.975"])),
         col = rgb(0, 0, 1, alpha = 0.2), border = NA)
} # }