1 #'@title First Differences of a Bayesian Logit or Probit model  2 #'@description R function to calculate first differences after a Bayesian logit or probit model.  3 #'First differences are a method to summarize effects across covariates. This quantity represents  4 #'the difference in predicted probabilities for each covariate for cases with low and high values  5 #'of the respective covariate. For each of these differences, all other variables are held constant  6 #'at their median. For more, see Long (1997, Sage Publications) and King, Tomz, and Wittenberg (2000,  7 #'American Journal of Political Science 44(2): 347-361).  8 #'@param modelmatrix model matrix, including intercept (if the intercept is among the  9 #'parameters estimated in the model). Create with model.matrix(formula, data).  10 #'Note: the order of columns in the model matrix must correspond to the order of columns  11 #'in the matrix of posterior draws in the \code{mcmcout} argument. See the \code{mcmcout}  12 #'argument for more.  13 #'@param mcmcout posterior distributions of all logit coefficients,  14 #'in matrix form. This can be created from rstan, MCMCpack, R2jags, etc. and transformed  15 #'into a matrix using the function as.mcmc() from the coda package for \code{jags} class  16 #'objects, as.matrix() from base R for \code{mcmc}, \code{mcmc.list}, \code{stanreg}, and  17 #'\code{stanfit} class objects, and \code{object$sims.matrix} for \code{bugs} class objects.  18 #'Note: the order of columns in this matrix must correspond to the order of columns  19 #'in the model matrix. One can do this by examining the posterior distribution matrix and sorting the  20 #'variables in the order of this matrix when creating the model matrix. A useful function for sorting  21 #'column names containing both characters and numbers as  22 #'you create the matrix of posterior distributions is \code{mixedsort()} from the gtools package.  23 #'@param link type of generalized linear model; a character vector set to \code{"logit"} (default)  24 #'or \code{"probit"}.  25 #'@param ci the bounds of the credible interval. Default is \code{c(0.025, 0.975)} for the 95\%  26 #'credible interval.  27 #'@param percentiles values of each predictor for which the difference in Pr(y = 1)  28 #'is to be calculated. Default is \code{c(0.25, 0.75)}, which will calculate the difference  29 #'between Pr(y = 1) for the 25th percentile and 75th percentile of the predictor. For binary  30 #'predictors, the function automatically calculates the difference between Pr(y = 1)  31 #'for x = 0 and x = 1.  32 #'@param fullsims logical indicator of whether full object (based on all MCMC draws  33 #'rather than their average) will be returned. Default is \code{FALSE}.  34 #'@references  35 #'\itemize{  36 #'\item King, Gary, Michael Tomz, and Jason Wittenberg. 2000. “Making the Most of Statistical  37 #'Analyses: Improving Interpretation and Presentation.” American Journal of Political Science  38 #'44 (2): 347–61. http://www.jstor.org/stable/2669316  39 #'\item Long, J. Scott. 1997. Regression Models for Categorical and Limited Dependent Variables.  40 #'Thousand Oaks: Sage Publications  41 #'}  42 #'@return An object of class \code{mcmcFD}. If \code{fullsims = FALSE} (default),  43 #'a data frame with five columns:  44 #'\itemize{  45 #'\item median_fd: median first difference  46 #'\item lower_fd: lower bound of credible interval of the first difference  47 #'\item upper_fd: upper bound of credible interval of the first difference  48 #'\item VarName: name of the variable as found in \code{modelmatrix}  49 #'\item VarID: identifier of the variable, based on the order of columns in  50 #'\code{modelmatrix} and \code{mcmcout}. Can be adjusted for plotting  51 #'}  52 #'If \code{fullsims = TRUE}, a matrix with as many columns as predictors in the model. Each row  53 #'is the first difference for that variable based on one set of posterior draws. Column names are taken  54 #'from the column names of \code{modelmatrix}.  55 #'  56 #'@examples  57 #' \dontshow{.old_wd <- setwd(tempdir())}  58 #' \donttest{  59 #' ## simulating data  60 #' set.seed(123456)  61 #' b0 <- 0.2 # true value for the intercept  62 #' b1 <- 0.5 # true value for first beta  63 #' b2 <- 0.7 # true value for second beta  64 #' n <- 500 # sample size  65 #' X1 <- runif(n, -1, 1)  66 #' X2 <- runif(n, -1, 1)  67 #' Z <- b0 + b1 * X1 + b2 * X2  68 #' pr <- 1 / (1 + exp(-Z)) # inv logit function  69 #' Y <- rbinom(n, 1, pr)  70 #' df <- data.frame(cbind(X1, X2, Y))  71 #'  72 #' ## formatting the data for jags  73 #' datjags <- as.list(df)  74 #' datjags$N <- length(datjags$Y)  75 #'  76 #' ## creating jags model  77 #' model <- function() {  78 #'  79 #' for(i in 1:N){  80 #' Y[i] ~ dbern(p[i]) ## Bernoulli distribution of y_i  81 #' logit(p[i]) <- mu[i] ## Logit link function  82 #' mu[i] <- b[1] +  83 #' b[2] * X1[i] +  84 #' b[3] * X2[i]  85 #' }  86 #'  87 #' for(j in 1:3){  88 #' b[j] ~ dnorm(0, 0.001) ## Use a coefficient vector for simplicity  89 #' }  90 #'  91 #' }  92 #'  93 #' params <- c("b")  94 #' inits1 <- list("b" = rep(0, 3))  95 #' inits2 <- list("b" = rep(0, 3))  96 #' inits <- list(inits1, inits2)  97 #'  98 #' ## fitting the model with R2jags  99 #' set.seed(123)  100 #' fit <- R2jags::jags(data = datjags, inits = inits,  101 #' parameters.to.save = params, n.chains = 2, n.iter = 2000,  102 #' n.burnin = 1000, model.file = model)  103 #'  104 #' ## running function with logit  105 #' xmat <- model.matrix(Y ~ X1 + X2, data = df)  106 #' mcmc <- coda::as.mcmc(fit)  107 #' mcmc_mat <- as.matrix(mcmc)[, 1:ncol(xmat)]  108 #' object <- mcmcFD(modelmatrix = xmat,  109 #' mcmcout = mcmc_mat)  110 #' object  111 #' }  112 #'  113 #' \dontshow{setwd(.old_wd)}  114 #'@export  115 #'  116 mcmcFD <- function(modelmatrix,  117  mcmcout,  118  link = "logit",  119  ci = c(0.025, 0.975),  120  percentiles = c(0.25, 0.75),  121  fullsims = FALSE) {  122   123 1  if(missing(modelmatrix) | missing(mcmcout)){  124 1  stop("Please enter required arguments.")  125  }  126 127 1  fdmat <- matrix(NA, ncol = 3, nrow = ncol(modelmatrix) - 1)  128 1  colnames(fdmat) <- c("median_fd", "lower_fd", "upper_fd")  129 1  rownames(fdmat) <- colnames(modelmatrix)[-1]  130 131 1  fdfull <- matrix(rep(NA),  132 1  ncol = ncol(modelmatrix) - 1,  133 1  nrow = nrow(mcmcout),  134 1  byrow = TRUE)  135 1  colnames(fdfull) <- colnames(modelmatrix)[-1]  136 137 1  for (i in 2:ncol(modelmatrix)){  138 139 1  X <- matrix(rep(apply(X = modelmatrix,  140 1  MARGIN = 2,  141 1  FUN = function(x) median(x)),  142 1  times = 2),  143 1  nrow = 2,  144 1  byrow = TRUE)  145 1  X[, i] <- ifelse(length(unique(modelmatrix[, i])) == 2 &  146 1  range(modelmatrix[, i]) == c(0, 1), c(0, 1),  147 1  quantile(modelmatrix[, i], probs = percentiles))  148 149  # X[, i] <- quantile(modelmatrix[, i], probs = percentiles)  150 151 1  if(link == "logit") {  152 1  pp <- plogis(t(X %*% t(mcmcout)))  153 1  } else if (link == "probit") {  154 1  pp <- pnorm(t(X %*% t(mcmcout)))  155  } else {  156 1  stop("Please enter valid link argument.")  157  }  158 159 160 1  fd <- pp[, 2] - pp[, 1]  161 162 1  fdmat[i-1, 1] <- quantile(fd, probs = c(0.5))  163 1  fdmat[i-1, 2] <- quantile(fd, probs = c(ci[1]))  164 1  fdmat[i-1, 3] <- quantile(fd, probs = c(ci[2]))  165 166 1  fdfull[, i-1] <- fd  167 168  }  169 170 1  fddat <- as.data.frame(fdmat)  171 1  fddat$VarName <- rownames(fdmat)  172 1  fddat$VarID <- row(fdmat)[, 1]  173   174 1  if (fullsims) {  175 1  return(structure(fdfull, fullsims = TRUE, class = c("mcmcFD", "matrix")))  176  } else {  177 1  return(structure(fddat, fullsims = FALSE, class = c("mcmcFD", "data.frame")))  178  }  179 180 }  181 182 #'@export  183 print.mcmcFD <- function(x, ...) {  184 1  if (attr(x, "fullsims")) {  185 1  print.table(x)  186  } else {  187 1  print.data.frame(x)  188  }  189 }  190 191 #'@title Plot Method for First Differences from MCMC output  192 #'@description The \code{plot} method for first differences generated from MCMC  193 #'output by \code{\link{mcmcFD}}. For more on this method, see Long  194 #'(1997, Sage Publications), and King, Tomz, and Wittenberg (2000, American  195 #'Journal of Political Science 44(2): 347-361). For a description of this type  196 #'of plot, see Figure 1 in Karreth (2018, International Interactions 44(3): 463-90).  197 #'@param x Output generated from \code{mcmcFD(..., full_sims = TRUE)}.  198 #'@param ROPE defaults to NULL. If not NULL, a numeric vector of length two,  199 #'defining the Region of Practical Equivalence around 0. See Kruschke (2013, Journal of  200 #'Experimental Psychology 143(2): 573-603) for more on the ROPE.  201 #'@param ... optional arguments to \code{\link[ggplot2]{theme}} from \code{\link[ggplot2:ggplot2-package]{ggplot2}}.  202 #'@references  203 #'\itemize{  204 #'\item Karreth, Johannes. 2018. “The Economic Leverage of International Organizations in Interstate Disputes.”  205 #'International Interactions 44 (3): 463-90. https://doi.org/10.1080/03050629.2018.1389728.  206 #'\item King, Gary, Michael Tomz, and Jason Wittenberg. 2000. “Making the Most of Statistical  207 #'Analyses: Improving Interpretation and Presentation.” American Journal of Political Science  208 #'44 (2): 347–61. http://www.jstor.org/stable/2669316.  209 #'\item Kruschke, John K. 2013. “Bayesian Estimation Supersedes the T-Test.” Journal of  210 #'Experimental Psychology: General 142 (2): 573–603. https://doi.org/10.1037/a0029146.  211 #'\item Long, J. Scott. 1997. Regression Models for Categorical and Limited Dependent Variables.  212 #'Thousand Oaks: Sage Publications.  213 #'}  214 #'@return a density plot of the differences in probabilities. The plot is made with ggplot2 and can be  215 #'passed on as an object to customize. Annotated numbers show the percent of posterior draws with  216 #'the same sign as the median estimate (if \code{ROPE = NULL}) or on the same side of the  217 #'ROPE as the median estimate (if \code{ROPE} is specified).  218 #'  219 #'@seealso \code{\link{mcmcFD}}  220 #'  221 #'@method plot mcmcFD  222 #'  223 #'@examples  224 #' \dontshow{.old_wd <- setwd(tempdir())}  225 #' \donttest{  226 #' ## simulating data  227 #' set.seed(123456)  228 #' b0 <- 0.2 # true value for the intercept  229 #' b1 <- 0.5 # true value for first beta  230 #' b2 <- 0.7 # true value for second beta  231 #' n <- 500 # sample size  232 #' X1 <- runif(n, -1, 1)  233 #' X2 <- runif(n, -1, 1)  234 #' Z <- b0 + b1 * X1 + b2 * X2  235 #' pr <- 1 / (1 + exp(-Z)) # inv logit function  236 #' Y <- rbinom(n, 1, pr)  237 #' df <- data.frame(cbind(X1, X2, Y))  238 #'  239 #' ## formatting the data for jags  240 #' datjags <- as.list(df)  241 #' datjags$N <- length(datjags$Y)  242 #'  243 #' ## creating jags model  244 #' model <- function() {  245 #'  246 #' for(i in 1:N){  247 #' Y[i] ~ dbern(p[i]) ## Bernoulli distribution of y_i  248 #' logit(p[i]) <- mu[i] ## Logit link function  249 #' mu[i] <- b[1] +  250 #' b[2] * X1[i] +  251 #' b[3] * X2[i]  252 #' }  253 #'  254 #' for(j in 1:3){  255 #' b[j] ~ dnorm(0, 0.001) ## Use a coefficient vector for simplicity  256 #' }  257 #'  258 #' }  259 #'  260 #' params <- c("b")  261 #' inits1 <- list("b" = rep(0, 3))  262 #' inits2 <- list("b" = rep(0, 3))  263 #' inits <- list(inits1, inits2)  264 #'  265 #' ## fitting the model with R2jags  266 #' set.seed(123)  267 #' fit <- R2jags::jags(data = datjags, inits = inits,  268 #' parameters.to.save = params, n.chains = 2, n.iter = 2000,  269 #' n.burnin = 1000, model.file = model)  270 #'  271 #' ## preparing data for mcmcFD()  272 #' xmat <- model.matrix(Y ~ X1 + X2, data = df)  273 #' mcmc <- coda::as.mcmc(fit)  274 #' mcmc_mat <- as.matrix(mcmc)[, 1:ncol(xmat)]  275 #'  276 #' ## plotting with mcmcFDplot()  277 #' full <- mcmcFD(modelmatrix = xmat,  278 #' mcmcout = mcmc_mat,  279 #' fullsims = TRUE)  280 #' plot(full)  281 #'  282 #' }  283 #'  284 #' \dontshow{setwd(.old_wd)}  285 #'  286 #' @export  287 #'  288 plot.mcmcFD <- function(x,  289  ROPE = NULL,  290  ...) {  291   292 1  if (!attr(x, "fullsims")) {  293 1  stop("full simulations must be used to plot posterior distribution")  294  }  295   296 1  ROPE <- check_ROPE_argument(ROPE)  297   298  ## multiply by 100 for percentage point change in output  299 1  x <- x * 100  300   301  # convert x to long data frame  302 1  fd_dat <- tidyr::gather(as.data.frame(x))  303   304  # create first plot  305   306 1  if(!is.null(ROPE)) {  307 1  fd_plot <- ggplot2::ggplot(data = fd_dat, aes(x = .data$value, y = .data$key)) +  308 1  ggplot2::geom_rect(xmin = ROPE[1], xmax = ROPE[2], ymin = 0, ymax = Inf, fill = "black") +  309 1  ggridges::stat_density_ridges(quantile_lines = TRUE,  310 1  quantiles = c(0.025, 0.5, 0.975),  311 1  vline_color = "white") +  312 1  ggplot2::scale_x_continuous(labels = function(x) x * 100) +  313 1  ggplot2::xlab("Percentage point change in Pr(y = 1)") +  314 1  ggplot2::ylab("")  315   316  # calculate area left/right of ROPE  317 1  fd_outROPE <- apply(x, 2,  318 1  function(x) ifelse(median(x) < 0,  319 1  sum(x < ROPE[1]) / length(x),  320 1  sum(x > ROPE[2]) / length(x)))  321 1  fd_annotate <- data.frame(xpos = apply(x, 2,  322 1  function(x) ifelse(median(x) < 0,  323 1  quantile(x, probs = 0.01) - 0.02,  324 1  quantile(x, probs = 0.99) + 0.02)),  325 1  ypos = as.factor(colnames(x)),  326 1  outROPE = paste(round(fd_outROPE * 100, digits = 1), "%", sep = ""))  327   328  # final plot  329 1  fd_plot <- fd_plot +  330 1  geom_text(data = fd_annotate, aes(x = .data$xpos, y = .data$ypos, label = .data$outROPE),  331 1  color = "black", nudge_y = 0.1, size = 4) +  332 1  ggplot2::theme(...)  333  } else {  334 1  fd_plot <- ggplot2::ggplot(data = fd_dat,  335 1  aes(x = .data$value, y = .data$key)) +  336 1  ggplot2::geom_vline(xintercept = 0) +  337 1  ggridges::stat_density_ridges(quantile_lines = TRUE,  338 1  quantiles = c(0.025, 0.5, 0.975),  339 1  vline_color = "white") +  340 1  ggplot2::xlab("Percentage point change in Pr(y = 1)") +  341 1  ggplot2::ylab("")  342   343  # calculate area left/right of 0  344 1  fd_out0 <- apply(x, 2,  345 1  function(x) ifelse(median(x) < 0,  346 1  sum(x < 0) / length(x),  347 1  sum(x > 0) / length(x)))  348 1  fd_annotate <- data.frame(xpos = apply(x, 2,  349 1  function(x) ifelse(median(x) < 0,  350 1  quantile(x, probs = 0.01) - 0.02,  351 1  quantile(x, probs = 0.99) + 0.02)),  352 1  ypos = as.factor(colnames(x)),  353 1  out0 = paste(round(fd_out0 * 100, digits = 1), "%", sep = ""))  354   355  # final plot  356 1  fd_plot <- fd_plot +  357 1  ggplot2::geom_text(data = fd_annotate,  358 1  aes(x = .data$xpos, y = .data$ypos, label = .data$out0),  359 1  color = "black", nudge_y = 0.1, size = 4) +  360 1  ggplot2::theme(...)  361   362  }  363   364 1  fd_plot  365   366 }  367 368 369 ## mcmcFDplot - deprecated  370 371 #'@title Plot First Differences from MCMC output  372 #'@description R function to plot first differences generated from MCMC output.  373 #'For more on this method, see the documentation for \code{mcmcFD()}, Long (1997,  374 #'Sage Publications), and King, Tomz, and Wittenberg (2000, American Journal  375 #'of Political Science 44(2): 347-361). For a description of this type of plot,  376 #'see Figure 1 in Karreth (2018, International Interactions 44(3): 463-90).  377 #'@param fdfull Output generated from \code{mcmcFD(..., full_sims = TRUE)}.  378 #'@param ROPE defaults to NULL. If not NULL, a numeric vector of length two,  379 #'defining the Region of Practical Equivalence around 0. See Kruschke (2013, Journal of  380 #'Experimental Psychology 143(2): 573-603) for more on the ROPE.  381 #'@references  382 #'\itemize{  383 #'\item Karreth, Johannes. 2018. “The Economic Leverage of International Organizations in Interstate Disputes.”  384 #'International Interactions 44 (3): 463-90. https://doi.org/10.1080/03050629.2018.1389728.  385 #'\item King, Gary, Michael Tomz, and Jason Wittenberg. 2000. “Making the Most of Statistical  386 #'Analyses: Improving Interpretation and Presentation.” American Journal of Political Science  387 #'44 (2): 347–61. http://www.jstor.org/stable/2669316.  388 #'\item Kruschke, John K. 2013. “Bayesian Estimation Supersedes the T-Test.” Journal of  389 #'Experimental Psychology: General 142 (2): 573–603. https://doi.org/10.1037/a0029146.  390 #'\item Long, J. Scott. 1997. Regression Models for Categorical and Limited Dependent Variables.  391 #'Thousand Oaks: Sage Publications.  392 #'}  393 #'@return a density plot of the differences in probabilities. The plot is made with ggplot2 and can be  394 #'passed on as an object to customize. Annotated numbers show the percent of posterior draws with  395 #'the same sign as the median estimate (if \code{ROPE = NULL}) or on the same side of the  396 #'ROPE as the median estimate (if \code{ROPE} is specified).  397 #'  398 #'@examples  399 #' \dontshow{.old_wd <- setwd(tempdir())}  400 #' \donttest{  401 #' ## simulating data  402 #' set.seed(123456)  403 #' b0 <- 0.2 # true value for the intercept  404 #' b1 <- 0.5 # true value for first beta  405 #' b2 <- 0.7 # true value for second beta  406 #' n <- 500 # sample size  407 #' X1 <- runif(n, -1, 1)  408 #' X2 <- runif(n, -1, 1)  409 #' Z <- b0 + b1 * X1 + b2 * X2  410 #' pr <- 1 / (1 + exp(-Z)) # inv logit function  411 #' Y <- rbinom(n, 1, pr)  412 #' df <- data.frame(cbind(X1, X2, Y))  413 #'  414 #' ## formatting the data for jags  415 #' datjags <- as.list(df)  416 #' datjags$N <- length(datjags$Y)  417 #'  418 #' ## creating jags model  419 #' model <- function() {  420 #'  421 #' for(i in 1:N){  422 #' Y[i] ~ dbern(p[i]) ## Bernoulli distribution of y_i  423 #' logit(p[i]) <- mu[i] ## Logit link function  424 #' mu[i] <- b[1] +  425 #' b[2] * X1[i] +  426 #' b[3] * X2[i]  427 #' }  428 #'  429 #' for(j in 1:3){  430 #' b[j] ~ dnorm(0, 0.001) ## Use a coefficient vector for simplicity  431 #' }  432 #'  433 #' }  434 #'  435 #' params <- c("b")  436 #' inits1 <- list("b" = rep(0, 3))  437 #' inits2 <- list("b" = rep(0, 3))  438 #' inits <- list(inits1, inits2)  439 #'  440 #' ## fitting the model with R2jags  441 #' set.seed(123)  442 #' fit <- R2jags::jags(data = datjags, inits = inits,  443 #' parameters.to.save = params, n.chains = 2, n.iter = 2000,  444 #' n.burnin = 1000, model.file = model)  445 #'  446 #' ## preparing data for mcmcFD()  447 #' xmat <- model.matrix(Y ~ X1 + X2, data = df)  448 #' mcmc <- coda::as.mcmc(fit)  449 #' mcmc_mat <- as.matrix(mcmc)[, 1:ncol(xmat)]  450 #'  451 #' ## plotting with mcmcFDplot()  452 #' full <- mcmcFD(modelmatrix = xmat,  453 #' mcmcout = mcmc_mat,  454 #' fullsims = TRUE)  455 #' # suppress deprecated warning for R check  456 #' suppressWarnings(mcmcFDplot(full))  457 #'  458 #' }  459 #'  460 #' \dontshow{setwd(.old_wd)}  461 #'  462 #' @name mcmcFDplot-deprecated  463 #' @usage mcmcFDplot(fdfull, ROPE = NULL)  464 #' @seealso \code{\link{BayesPostEst-deprecated}}  465 #' @seealso \code{\link{mcmcFD}}  466 #' @keywords internal  467 #'  468 #' @rdname BayesPostEst-deprecated  469 #' @section \code{mcmcFDplot}:  470 #' For \code{mcmcFDplot}, use \code{\link{plot.mcmcFD}}.  471 #'  472 #' @export  473 #'  474 475 mcmcFDplot <- function(fdfull,  476  ROPE = NULL){  477   478 1  .Deprecated("plot.mcmcFD", package = "BayesPostEst")  479   480 1  ROPE <- check_ROPE_argument(ROPE)  481   482 1  fdfull <- fdfull * 100  483   484 1  fd_dat <- tidyr::gather(as.data.frame(fdfull))  485   486  # create first plot  487   488 1  if(!is.null(ROPE)) {  489 1  fd_plot <- ggplot2::ggplot(data = fd_dat, aes(x = .data$value, y = .data$key)) +  490 1  ggplot2::geom_rect(xmin = ROPE[1], xmax = ROPE[2], ymin = 0, ymax = Inf, fill = "black") +  491 1  ggridges::stat_density_ridges(quantile_lines = TRUE,  492 1  quantiles = c(0.025, 0.5, 0.975),  493 1  vline_color = "white") +  494 1  ggplot2::scale_x_continuous(labels = function(x) x * 100) +  495 1  ggplot2::xlab("Percentage point change in Pr(y = 1)") +  496 1  ggplot2::ylab("")  497   498  # calculate area left/right of ROPE  499 1  fd_outROPE <- apply(fdfull, 2,  500 1  function(x) ifelse(median(x) < 0,  501 1  sum(x < ROPE[1]) / length(x),  502 1  sum(x > ROPE[2]) / length(x)))  503 1  fd_annotate <- data.frame(xpos = apply(fdfull, 2,  504 1  function(x) ifelse(median(x) < 0,  505 1  quantile(x, probs = 0.01) - 0.02,  506 1  quantile(x, probs = 0.99) + 0.02)),  507 1  ypos = as.factor(colnames(fdfull)),  508 1  outROPE = paste(round(fd_outROPE * 100, digits = 1), "%", sep = ""))  509   510  # final plot  511 1  fd_plot <- fd_plot +  512 1  geom_text(data = fd_annotate, aes(x = .data$xpos, y = .data$ypos, label = .data$outROPE),  513 1  color = "black", nudge_y = 0.1, size = 4)  514  } else {  515 1  fd_plot <- ggplot2::ggplot(data = fd_dat,  516 1  aes(x = .data$value, y = .data$key)) +  517 1  ggplot2::geom_vline(xintercept = 0) +  518 1  ggridges::stat_density_ridges(quantile_lines = TRUE,  519 1  quantiles = c(0.025, 0.5, 0.975),  520 1  vline_color = "white") +  521 1  ggplot2::scale_x_continuous(labels = function(x) x*100) +  522 1  ggplot2::xlab("Percentage point change in Pr(y = 1)") +  523 1  ggplot2::ylab("")  524   525  # calculate area left/right of 0  526 1  fd_out0 <- apply(fdfull, 2,  527 1  function(x) ifelse(median(x) < 0,  528 1  sum(x < 0) / length(x),  529 1  sum(x > 0) / length(x)))  530 1  fd_annotate <- data.frame(xpos = apply(fdfull, 2,  531 1  function(x) ifelse(median(x) < 0,  532 1  quantile(x, probs = 0.01) - 0.02,  533 1  quantile(x, probs = 0.99) + 0.02)),  534 1  ypos = as.factor(colnames(fdfull)),  535 1  out0 = paste(round(fd_out0 * 100, digits = 1), "%", sep = ""))  536   537  # final plot  538 1  fd_plot <- fd_plot +  539 1  ggplot2::geom_text(data = fd_annotate,  540 1  aes(x = .data$xpos, y = .data$ypos, label = .data\$out0),  541 1  color = "black", nudge_y = 0.1, size = 4)  542  }  543   544 1  fd_plot  545   546 } 

Read our documentation on viewing source code .