`R/mcmcAveProb.R`

`mcmcAveProb.Rd`

This function calculates predicted probabilities for "average" cases after a Bayesian logit or probit model. As "average" cases, this function calculates the median value of each predictor. For an explanation of predicted probabilities for "average" cases, see e.g. King, Tomz & Wittenberg (2000, American Journal of Political Science 44(2): 347-361).

mcmcAveProb(modelmatrix, mcmcout, xcol, xrange, xinterest, link = "logit", ci = c(0.025, 0.975), fullsims = FALSE)

modelmatrix | model matrix, including intercept (if the intercept is among the
parameters estimated in the model). Create with model.matrix(formula, data).
Note: the order of columns in the model matrix must correspond to the order of columns
in the matrix of posterior draws in the |
---|---|

mcmcout | posterior distributions of all logit coefficients,
in matrix form. This can be created from rstan, MCMCpack, R2jags, etc. and transformed
into a matrix using the function as.mcmc() from the coda package for |

xcol | column number of the posterior draws ( |

xrange | name of the vector with the range of relevant values of the explanatory variable for which to calculate associated Pr(y = 1). |

xinterest | semi-optional argument. Name of the explanatory variable for which
to calculate associated Pr(y = 1). If |

link | type of generalized linear model; a character vector set to |

ci | the bounds of the credible interval. Default is |

fullsims | logical indicator of whether full object (based on all MCMC draws
rather than their average) will be returned. Default is |

if `fullsims = FALSE`

(default), a tibble with 4 columns:

x: value of variable of interest, drawn from

`xrange`

median_pp: median predicted Pr(y = 1) when variable of interest is set to x, holding all other predictors to average (median) values

lower_pp: lower bound of credible interval of predicted probability at given x

upper_pp: upper bound of credible interval of predicted probability at given x

if `fullsims = TRUE`

, a tibble with 3 columns:

Iteration: number of the posterior draw

x: value of variable of interest, drawn from

`xrange`

pp: average predicted Pr(y = 1) when variable of interest is set to x, holding all other predictors to average (median) values

This function calculates predicted probabilities for "average" cases after a Bayesian logit or probit model. For an explanation of predicted probabilities for "average" cases, see e.g. King, Tomz & Wittenberg (2000, American Journal of Political Science 44(2): 347-361)

King, Gary, Michael Tomz, and Jason Wittenberg. 2000. “Making the Most of Statistical Analyses: Improving Interpretation and Presentation.” American Journal of Political Science 44 (2): 347–61. http://www.jstor.org/stable/2669316

.old_wd <- setwd(tempdir()) # \donttest{ ## 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 pr <- 1 / (1 + exp(-Z)) # inv logit function Y <- rbinom(n, 1, pr) data <- data.frame(cbind(X1, X2, Y)) ## formatting the data for jags datjags <- as.list(data) datjags$N <- length(datjags$Y) ## creating jags model model <- function() { for(i in 1:N){ Y[i] ~ dbern(p[i]) ## Bernoulli distribution of y_i logit(p[i]) <- mu[i] ## Logit link function mu[i] <- b[1] + b[2] * X1[i] + b[3] * X2[i] } for(j in 1:3){ b[j] ~ dnorm(0, 0.001) ## Use a coefficient vector for simplicity } } params <- c("b") inits1 <- list("b" = rep(0, 3)) inits2 <- list("b" = rep(0, 3)) inits <- list(inits1, inits2) ## fitting the model with R2jags library(R2jags)#>#>#>#>#> #>#>#> #>set.seed(123) fit <- jags(data = datjags, inits = inits, parameters.to.save = params, n.chains = 2, n.iter = 2000, n.burnin = 1000, model.file = model)#> Compiling model graph #> Resolving undeclared variables #> Allocating nodes #> Graph information: #> Observed stochastic nodes: 500 #> Unobserved stochastic nodes: 3 #> Total graph size: 3506 #> #> Initializing model #>### average value approach library(coda) xmat <- model.matrix(Y ~ X1 + X2, data = data) mcmc <- as.mcmc(fit) mcmc_mat <- as.matrix(mcmc)[, 1:ncol(xmat)] X1_sim <- seq(from = min(datjags$X1), to = max(datjags$X1), length.out = 10) ave_prob <- mcmcAveProb(modelmatrix = xmat, mcmcout = mcmc_mat, xrange = X1_sim, xcol = 2) # } setwd(.old_wd)