hmclearn (version 0.0.5)

predict.hmclearn: Model Predictions for HMC or MH

Description

predict generates simulated data from the posterior predictive distribution. This simulated data can be used for posterior predictive check diagnostics from the bayesplot package

Usage

# S3 method for hmclearn
predict(object, X, fam = "linear", burnin = NULL, draws = NULL, ...)

Arguments

object

an object of class hmclearn, usually a result of a call to mh or hmc

X

design matrix, either from fitting the model or new data

fam

generalized linear model family. Currently "linear", "binomial", and "poisson" are supported

burnin

optional numeric parameter for the number of initial MCMC samples to omit from the summary

draws

Number of simulated values from the posterior conditioned on X

...

additional parameters, currently unsupported

Value

An object of class hmclearnpred.

Elements of <code>hmclearnpred</code> objects

y

Median simulated values for each observation in X

yrep

Matrix of simulated values where each row is a draw from the posterior predictive distribution

X

Numeric design matrix

References

Gabry, Jonah and Mahr, Tristan (2019). bayesplot: Plotting for Bayesian Models. https://mc-stan.org/bayesplot/

Gabry, J., Simpson, D., Vehtari, A., Betancourt, M., and Gelman, A (2019). Visualization in Bayesian Workflow. Journal of the Royal Statistical Society: Series A. Vol 182. Issue 2. p.389-402.

Examples

# NOT RUN {
# Linear regression example
set.seed(521)
X <- cbind(1, matrix(rnorm(300), ncol=3))
betavals <- c(0.5, -1, 2, -3)
y <- X%*%betavals + rnorm(100, sd=.2)

f1 <- hmc(N = 500,
          theta.init = c(rep(0, 4), 1),
          epsilon = 0.01,
          L = 10,
          logPOSTERIOR = linear_posterior,
          glogPOSTERIOR = g_linear_posterior,
          varnames = c(paste0("beta", 0:3), "log_sigma_sq"),
          param=list(y=y, X=X), parallel=FALSE, chains=1)

summary(f1)

p <- predict(f1, X)
predvals <- p$y
plot(predvals, y, xlab="predicted", ylab="actual")

X2 <- cbind(1, matrix(rnorm(30), ncol=3))
p2 <- predict(f1, X2)
p2$y

# }