1
#' Marginal Effects Plots for MCMC Output
2
#'
3
#' Marginal effects plots for MCMC output using \code{ggplot2}
4
#'
5
#' @param mod Bayesian model object generated by R2jags, rjags, R2WinBUGS, R2OpenBUGS, 
6
#' MCMCpack, rstan, rstanarm, and brms.
7
#' @param main a character with the name of the parameter of interest in the
8
#' interaction term.
9
#' @param int a character with the name of the moderating parameter in the
10
#' interaction term.
11
#' @param moderator a vector of values that the moderating parameter takes on
12
#' in the data.
13
#' @param pointest a character indicating whether to use the mean or median for
14
#' point estimates in the plot.
15
#' @param seq a numeric giving the number of moderator values used to generate
16
#' the marginal effects plot.
17
#' @param ci a scalar indicating the confidence level of the uncertainty intervals.
18
#' @param hpdi a logical indicating whether to use highest posterior density intervals
19
#' or equal tailed credible intervals to capture uncertainty.
20
#' @param plot logical indicating whether to return a \code{ggplot} object or the
21
#' underlying tidy DataFrame. By default, \code{mcmcMargEff} returns a line and
22
#' ribbon plot for continuous variables, and a dot and line plot for factor
23
#' variables and discrete variables with fewer than 25 unique values.
24
#' @param xlab character giving x axis label if \code{plot = TRUE}, default \code{"Moderator"}
25
#' @param ylab character giving y axis label if \code{plot = TRUE}, default \code{"Marginal Effect"}
26
#'
27
#' @return a \code{ggplot} object or a tidy DataFrame.
28
#' 
29
#' @author Rob Williams, \email{jayrobwilliams@gmail.com}
30
#'
31
#' @examples
32
#' \dontshow{.old_wd <- setwd(tempdir())}
33
#' \donttest{
34
#' ## simulating data
35
#' set.seed(123456)
36
#' b0 <- 0.2 # true value for the intercept
37
#' b1 <- 0.5 # true value for first beta
38
#' b2 <- 0.7 # true value for second beta
39
#' n <- 500 # sample size
40
#' X1 <- runif(n, -1, 1)
41
#' X2 <- runif(n, -1, 1)
42
#' Z <- b0 + b1 * X1 + b2 * X2
43
#'
44
#' ## linear model data
45
#' Y_linear <- rnorm(n, Z, 1)
46
#' df <- data.frame(cbind(X1, X2, Y = Y_linear))
47
#' 
48
#' ## formatting the data for jags
49
#' datjags <- as.list(df)
50
#' datjags$N <- length(datjags$Y)
51
#' 
52
#' ## creating jags model
53
#' model <- function()  {
54
#'   
55
#'   for(i in 1:N){
56
#'     Y[i] ~ dnorm(mu[i], sigma)  ## Bernoulli distribution of y_i
57
#'     
58
#'     mu[i] <- b[1] + 
59
#'       b[2] * X1[i] + 
60
#'       b[3] * X2[i] +
61
#'       b[4] * X1[i] * X2[i]
62
#'     
63
#'   }
64
#'   
65
#'   for(j in 1:4){
66
#'     b[j] ~ dnorm(0, 0.001) ## Use a coefficient vector for simplicity
67
#'   }
68
#'   
69
#'   sigma ~ dexp(1)
70
#'   
71
#' }
72
#' 
73
#' params <- c("b")
74
#' inits1 <- list("b" = rep(0, 4))
75
#' inits2 <- list("b" = rep(0, 4))
76
#' inits <- list(inits1, inits2)
77
#' 
78
#' ## fitting the model with R2jags
79
#' set.seed(123)
80
#' fit <- R2jags::jags(data = datjags, inits = inits, 
81
#'                     parameters.to.save = params, n.chains = 2, n.iter = 2000, 
82
#'                     n.burnin = 1000, model.file = model)
83
#'
84
#' mcmcMargEff(mod = fit,
85
#'             main = 'b[2]',
86
#'             int = 'b[4]',
87
#'             moderator = sim_data_interactive$X2,
88
#'             plot = TRUE)
89
#' }
90
#' 
91
#' \dontshow{setwd(.old_wd)}
92
#' @export
93
#' 
94
mcmcMargEff <- function(mod, 
95
                        main, 
96
                        int, 
97
                        moderator, 
98
                        pointest = 'mean', 
99
                        seq = 100,
100
                        ci = .95, 
101
                        hpdi = FALSE, 
102
                        plot = TRUE, 
103
                        xlab = 'Moderator',
104
                        ylab = 'Marginal Effect') {
105

106
  ## pull in unexported functions from other packages
107
  ## other options for future versions might include lifting this and adding authors as copr holders
108 1
  runjags.as.mcmc.list.runjags = getFromNamespace("as.mcmc.list.runjags", "runjags")
109 1
  coda.hpdinterval.mcmc = getFromNamespace("HPDinterval.mcmc", "coda")
110
  
111 1
  if (inherits(mod, what = c("jags", "rjags"))) {
112 1
    samps <- as.matrix(coda::as.mcmc(mod))
113
  }
114 1
  if (inherits(mod, what = "bugs")) {
115 0
    samps <- mod$sims.matrix
116
  }
117 1
  if (inherits(mod, what = "runjags")) {
118 0
    samps <- as.matrix(runjags.as.mcmc.list.runjags(mod))
119
  }
120 1
  if (inherits(mod, what = c("mcmc", "mcmc.list", "stanfit", "stanreg",
121 1
                             "brmsfit"))) {
122 0
    samps <- as.matrix(mod)
123
  }
124
  
125 1
  samps <- samps[, c(main, int)]
126

127
  ## expand moderating variable to range of values
128 1
  if(!is.factor(moderator) & all(unique(moderator) %% 1 != 0)) {
129 1
    mod_range <- seq(min(moderator), max(moderator), length.out = seq)
130 1
    categorical <- F
131 1
  } else if ((is.factor(moderator) | all(unique(moderator) %% 1 == 0)) &
132 1
             length(unique(moderator)) >= 25) {
133 0
    mod_range <- seq(min(moderator), max(moderator), length.out = seq)
134 0
    categorical <- F
135 1
  } else if (is.factor(moderator) | all(unique(moderator) %% 1 == 0)) {
136 1
    mod_range <- sort(unique(moderator))
137 1
    seq = length(mod_range)
138 1
    categorical <- T
139
  }
140
  
141
  ## compute marginal effect for each sample
142 1
  marg <- rep(samps[, 1], seq) + samps[, 2] %o% mod_range
143
  
144 1
  if (pointest == 'mean') {
145 1
    marg_pe <- apply(marg, 2, mean)
146 1
  } else if (pointest == 'median') {
147 1
    marg_pe <- apply(marg, 2, median)
148
  } else {
149 0
    stop("pointest must be either 'mean' or 'median'")
150
  }
151
  
152
  ## calculate marginal effect for mean
153 1
  if (!hpdi) {
154 1
    marg_ci<- t(apply(marg, 2, quantile, probs = c(.5 - ci/2, .5 + ci/2)))
155 1
  } else if (hpdi) {
156 1
    marg_ci <- t(apply(marg, 2, coda.hpdinterval.mcmc, prob = ci))
157
  } else {
158 0
    stop("hpdi must be either true or false")
159
  }
160

161
  ## create dataframe for plotting
162 1
  marg_gg <- data.frame(mod = mod_range, pe = marg_pe,
163 1
                        lo = marg_ci[, 1], hi = marg_ci[, 2])
164
  
165
  ## return marginal effects plot or underlying dataframe
166 1
  if (!plot) {
167 1
    marg_gg
168
  } else {
169 1
    if (categorical) {
170 0
      ggplot2::ggplot(data = marg_gg,
171 0
                      ggplot2::aes(x = .data$mod, y = .data$pe,
172 0
                                   ymin = .data$lo, ymax = .data$hi)) +
173 0
        ggplot2::geom_linerange() +
174 0
        ggplot2::geom_hline(yintercept = 0, lty = 2, color = 'gray40', lwd = .5) +
175 0
        ggplot2::geom_point() +
176 0
        ggplot2::labs(x = xlab, y = ylab)
177
    } else {
178 1
      ggplot2::ggplot(data = marg_gg,
179 1
                      ggplot2::aes(x = .data$mod, y = .data$pe,
180 1
                                   ymin = .data$lo, ymax = .data$hi)) +
181 1
        ggplot2::geom_ribbon(alpha = .25) +
182 1
        ggplot2::geom_hline(yintercept = 0, lty = 2, color = 'gray40', lwd = .5) +
183 1
        ggplot2::geom_line() +
184 1
        ggplot2::labs(x = xlab, y = ylab)
185
    }
186
  }
187
}

Read our documentation on viewing source code .

Loading