1
#'This function calculates predicted probabilities for 
2
#'"observed" cases after a Bayesian logit or probit model
3
#'following Hanmer and Kalkan (2013, American Journal of 
4
#'Political Science 57(1): 263-277)
5
#'@title Predicted Probabilities using Bayesian MCMC estimates for the Average of Observed Cases
6
#'@description Implements R function to calculate the predicted probabilities
7
#'for "observed" cases after a Bayesian logit or probit model, following
8
#'Hanmer & Kalkan (2013) (2013, American Journal of Political Science 57(1): 263-277).
9
#'@param modelmatrix model matrix, including intercept (if the intercept is among the
10
#'parameters estimated in the model). Create with model.matrix(formula, data).
11
#'Note: the order of columns in the model matrix must correspond to the order of columns 
12
#'in the matrix of posterior draws in the \code{mcmcout} argument. See the \code{mcmcout}
13
#'argument for more.
14
#'@param mcmcout posterior distributions of all logit coefficients, 
15
#'in matrix form. This can be created from rstan, MCMCpack, R2jags, etc. and transformed
16
#'into a matrix using the function as.mcmc() from the coda package for \code{jags} class
17
#'objects, as.matrix() from base R for \code{mcmc}, \code{mcmc.list}, \code{stanreg}, and 
18
#'\code{stanfit} class objects, and \code{object$sims.matrix} for \code{bugs} class objects.
19
#'Note: the order of columns in this matrix must correspond to the order of columns 
20
#'in the model matrix. One can do this by examining the posterior distribution matrix and sorting the 
21
#'variables in the order of this matrix when creating the model matrix. A useful function for sorting 
22
#'column names containing both characters and numbers as 
23
#'you create the matrix of posterior distributions is \code{mixedsort()} from the gtools package.
24
#'@param xcol column number of the posterior draws (\code{mcmcout}) and model matrices 
25
#'that corresponds to the explanatory variable for which to calculate associated Pr(y = 1).
26
#'Note that the columns in these matrices must match.
27
#'@param xrange name of the vector with the range of relevant values of the 
28
#'explanatory variable for which to calculate associated Pr(y = 1).
29
#'@param xinterest semi-optional argument. Name of the explanatory variable for which 
30
#'to calculate associated Pr(y = 1). If \code{xcol} is supplied, this is not needed. 
31
#'If both are supplied, the function defaults to \code{xcol} and this argument is ignored.
32
#'@param link type of generalized linear model; a character vector set to \code{"logit"} (default) 
33
#'or \code{"probit"}.
34
#'@param ci the bounds of the credible interval. Default is \code{c(0.025, 0.975)} for the 95\% 
35
#'credible interval.
36
#'@param fullsims logical indicator of whether full object (based on all MCMC draws 
37
#'rather than their average) will be returned. Default is \code{FALSE}. Note: The longer 
38
#'\code{xrange} is, the larger the full output will be if \code{TRUE} is selected.
39
#'@references Hanmer, Michael J., & Ozan Kalkan, K. (2013). Behind the curve: Clarifying 
40
#'the best approach to calculating predicted probabilities and marginal effects from 
41
#'limited dependent variable models. American Journal of Political Science, 57(1), 
42
#'263-277. https://doi.org/10.1111/j.1540-5907.2012.00602.x
43
#'@return if \code{fullsims = FALSE} (default), a tibble with 4 columns:
44
#'\itemize{
45
#'\item x: value of variable of interest, drawn from \code{xrange}
46
#'\item median_pp: median predicted Pr(y = 1) when variable of interest is set to x
47
#'\item lower_pp: lower bound of credible interval of predicted probability at given x
48
#'\item upper_pp: upper bound of credible interval of predicted probability at given x
49
#'}
50
#'if \code{fullsims = TRUE}, a tibble with 3 columns:
51
#'\itemize{
52
#'\item Iteration: number of the posterior draw
53
#'\item x: value of variable of interest, drawn from \code{xrange}
54
#'\item pp: average predicted Pr(y = 1) of all observed cases when variable of interest is set to x
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
#' library(R2jags)
100
#' set.seed(123)
101
#' fit <- jags(data = datjags, inits = inits, 
102
#'           parameters.to.save = params, n.chains = 2, n.iter = 2000, 
103
#'           n.burnin = 1000, model.file = model)
104
#' 
105
#' ### observed value approach
106
#' library(coda)
107
#' xmat <- model.matrix(Y ~ X1 + X2, data = df)
108
#' mcmc <- as.mcmc(fit)
109
#' mcmc_mat <- as.matrix(mcmc)[, 1:ncol(xmat)]
110
#' X1_sim <- seq(from = min(datjags$X1),
111
#'               to = max(datjags$X1), 
112
#'               length.out = 10)
113
#' obs_prob <- mcmcObsProb(modelmatrix = xmat,
114
#'                         mcmcout = mcmc_mat,
115
#'                         xrange = X1_sim,
116
#'                         xcol = 2)
117
#' }
118
#' 
119
#' \dontshow{setwd(.old_wd)}
120
#'@export
121
#'
122
mcmcObsProb <- function(modelmatrix,
123
                        mcmcout, 
124
                        xcol, 
125
                        xrange, 
126
                        xinterest,
127
                        link = "logit", 
128
                        ci = c(0.025, 0.975),
129
                        fullsims = FALSE){
130
  
131
  # checking arguments
132 1
  if(missing(xcol) & missing(xinterest)) {
133 0
    stop("Please enter a column number or name of your variable of interest)")
134
  }
135 1
  if(!missing(xcol) & !missing(xinterest)) {
136 0
    message("Both xcol and xinterest were supplied by user. Function defaults to xcol")
137
  }
138 1
  if(!missing(xinterest)) {
139 0
    if(!(xinterest %in% variable.names(modelmatrix)))
140 0
      stop("Variable name does not match any in the matrix. Please enter another.")
141
  }
142
  
143 1
  X <- matrix(rep(t(modelmatrix), length(xrange)), 
144 1
              ncol = ncol(modelmatrix), byrow = TRUE )
145 1
  colnames(X) <- variable.names(modelmatrix)
146 1
  if(!missing(xcol)) {
147 1
    X[, xcol] <- sort(rep(xrange, times = nrow(X) / length(xrange)))
148
  } else {
149 0
    X[ , grepl( xinterest , variable.names( X ) ) ] <- 
150 0
      sort(rep(xrange, times = nrow(X) / length(xrange)))
151
  }
152
  
153 1
  if(link == "logit"){
154 1
    pp <- plogis(t(X %*% t(mcmcout)))
155
  }
156
  
157 1
  if(link == "probit"){
158 1
    pp <- pnorm(t(X %*% t(mcmcout)))
159
  }
160
  
161
  
162
  # emptry matrix for PPs
163 1
  pp_mat <- matrix(NA, nrow = nrow(mcmcout), ncol = length(xrange))
164
  
165
  # indices
166 1
  pp_mat_lowerindex <- 1 + (0:(length(xrange) - 1) * nrow(modelmatrix))
167 1
  pp_mat_upperindex <- nrow(modelmatrix) + (0:(length(xrange) - 1) * 
168 1
                                               nrow(modelmatrix))
169
  
170
  
171
  # fill matrix with PPs, one for each value of the predictor of interest
172 1
  for(i in 1:length(xrange)){
173 1
    pp_mat[, i] <- apply(X = pp[, 
174 1
                                c(pp_mat_lowerindex[i]:pp_mat_upperindex[i])], 
175 1
                         MARGIN = 1, FUN = function(x) mean(x))
176
  }
177
  
178 1
  median_pp <- apply(X = pp_mat, MARGIN = 2, function(x) quantile(x, probs = c(0.5)))
179 1
  lower_pp <- apply(X = pp_mat, MARGIN = 2, function(x) quantile(x, probs = ci[1]))
180 1
  upper_pp <- apply(X = pp_mat, MARGIN = 2, function(x) quantile(x, probs = ci[2]))
181
  
182 1
  pp_dat <- dplyr::tibble(x = xrange,
183 1
                   median_pp = median_pp,
184 1
                   lower_pp = lower_pp,
185 1
                   upper_pp = upper_pp)
186
  
187 1
  if(fullsims == FALSE){
188 1
    return(pp_dat) # pp_dat was created by summarizing longFrame
189
  }
190
  
191 1
  if(fullsims == TRUE){
192 1
    longFrame <- reshape2::melt(pp_mat, id.vars = .data$Var2)
193 1
    names(longFrame) <- c("Iteration", "x", "pp")
194 1
    return(longFrame) 
195
  }
196
  
197
}

Read our documentation on viewing source code .

Loading