Marginal Effects Plots for MCMC Output

Marginal effects plots for MCMC output using ggplot2

mcmcMargEff(
  mod,
  main,
  int,
  moderator,
  pointest = "mean",
  seq = 100,
  ci = 0.95,
  hpdi = FALSE,
  plot = TRUE,
  xlab = "Moderator",
  ylab = "Marginal Effect"
)

Arguments

mod	Bayesian model object generated by R2jags, rjags, R2WinBUGS, R2OpenBUGS, MCMCpack, rstan, rstanarm, and brms.
main	a character with the name of the parameter of interest in the interaction term.
int	a character with the name of the moderating parameter in the interaction term.
moderator	a vector of values that the moderating parameter takes on in the data.
pointest	a character indicating whether to use the mean or median for point estimates in the plot.
seq	a numeric giving the number of moderator values used to generate the marginal effects plot.
ci	a scalar indicating the confidence level of the uncertainty intervals.
hpdi	a logical indicating whether to use highest posterior density intervals or equal tailed credible intervals to capture uncertainty.
plot	logical indicating whether to return a ggplot object or the underlying tidy DataFrame. By default, mcmcMargEff returns a line and ribbon plot for continuous variables, and a dot and line plot for factor variables and discrete variables with fewer than 25 unique values.
xlab	character giving x axis label if plot = TRUE, default "Moderator"
ylab	character giving y axis label if plot = TRUE, default "Marginal Effect"

Value

a ggplot object or a tidy DataFrame.

Author

Rob Williams, jayrobwilliams@gmail.com

Examples

.old_wd <- setwd(tempdir())
# \donttest{
if (interactive()) {
## simulating data
set.seed(123456)
b0 <- 0.2 # true value for the intercept
b1 <- 0.5 # true value for first beta
b2 <- 0.7 # true value for second beta
n <- 500 # sample size
X1 <- runif(n, -1, 1)
X2 <- runif(n, -1, 1)
Z <- b0 + b1 * X1 + b2 * X2

## linear model data
Y_linear <- rnorm(n, Z, 1)
df <- data.frame(cbind(X1, X2, Y = Y_linear))

## formatting the data for jags
datjags <- as.list(df)
datjags$N <- length(datjags$Y)

## creating jags model
model <- function()  {
  
  for(i in 1:N){
    Y[i] ~ dnorm(mu[i], sigma)  ## Bernoulli distribution of y_i
    
    mu[i] <- b[1] + 
      b[2] * X1[i] + 
      b[3] * X2[i] +
      b[4] * X1[i] * X2[i]
    
  }
  
  for(j in 1:4){
    b[j] ~ dnorm(0, 0.001) ## Use a coefficient vector for simplicity
  }
  
  sigma ~ dexp(1)
  
}

params <- c("b")
inits1 <- list("b" = rep(0, 4))
inits2 <- list("b" = rep(0, 4))
inits <- list(inits1, inits2)

## fitting the model with R2jags
set.seed(123)
fit <- R2jags::jags(data = datjags, inits = inits, 
                    parameters.to.save = params, n.chains = 2, n.iter = 2000, 
                    n.burnin = 1000, model.file = model)

mcmcMargEff(mod = fit,
            main = 'b[2]',
            int = 'b[4]',
            moderator = sim_data_interactive$X2,
            plot = TRUE)
}
# }

setwd(.old_wd)