1
#' @title LaTeX or HTML regression tables for MCMC Output
2
#' @description This function creates LaTeX or HTML regression tables for MCMC Output using 
3
#' the \code{\link[texreg]{texreg}} function from the \code{\link[texreg:texreg-package]{texreg}} R package.
4
#' @param mod Bayesian model object generated by R2jags, rjags, R2WinBUGS, R2OpenBUGS, 
5
#' MCMCpack, rstan, rstanarm, and brms, or a list of model objects of the same class.
6
#' @param pars a scalar or vector of the parameters you wish to include in the table.
7
#' By default, \code{mcmcReg} includes all parameters saved in a model object. If a
8
#' model has lots of samples and lots of saved parameters, not explicitly specifying
9
#' a limited number of parameters to include via \code{pars} may take a long time.
10
#' \code{pars} can either be a vector with the specific parameters to be included
11
#' in the table e.g. \code{pars = c("beta[1]", "beta[2]", "beta[3]")}, or they can
12
#' be partial names that will be matched using regular expressions e.g.
13
#' \code{pars = "beta"} if \code{regex = TRUE}. Both of these will include
14
#' \code{beta[1]}, \code{beta[2]}, and \code{beta[3]} in the table. When
15
#' combining models with different parameters in one table, this argument also
16
#' accepts a list the length of the number of models.
17
#' @param pointest a character indicating whether to use the mean or median for
18
#' point estimates in the table.
19
#' @param ci a scalar indicating the confidence level of the uncertainty intervals.
20
#' @param hpdi a logical indicating whether to use highest posterior density
21
#' intervals instead of equal tailed credible intervals to capture uncertainty
22
#' (default \code{FALSE}).
23
#' @param sd a logical indicating whether to report the standard deviation of
24
#' posterior distributions instead of an uncertainty interval
25
#' (default \code{FALSE}). If \code{TRUE}, overrides \code{ci}, \code{hpdi}, and
26
#' \code{pr}.
27
#' @param pr a logical indicating whether to report the probability that a
28
#' coefficient is in the same direction as the point estimate for that
29
#' coefficient (default \code{FALSE}). If \code{TRUE}, overrides \code{ci} and
30
#' \code{hpdi}.
31
#' @param coefnames an optional vector or list of vectors containing parameter
32
#' names for each model. If there are multiple models, the list must have the same
33
#' number of elements as there are models, and the vector of names in each list
34
#' element must match the number of parameters. If not supplied, the function
35
#' will use the parameter names in the model object(s). Note that this replaces
36
#' the standard \code{custom.coef.names} argument in \code{\link[texreg]{texreg}}
37
#' because there is no \code{extract} method for MCMC model objects, and many
38
#' MCMC model objects do not have unique parameter names.
39
#' @param gof a named list of goodness of fit statistics, or a list of such lists.
40
#' @param gofnames an optional vector or list of vectors containing
41
#' goodness of fit statistic names for each model. Like \code{coefnames} in this function
42
#' (which replaces the \code{custom.coef.names} argument in \code{\link[texreg]{texreg}}),
43
#' \code{gofnames} replaces the standard \code{custom.gof.names} argument in
44
#' \code{\link[texreg]{texreg}}. If 
45
#' there are multiple models, the list must have the same number of elements as
46
#' there are models, and the vector of names in each list element must match the
47
#' number of goodness of fit statistics.
48
#' @param format a character indicating \code{latex} or \code{html} output.
49
#' @param file optional file name to write table to file instead of printing to
50
#' console.
51
#' @param regex use regular expression matching with \code{pars}?
52
#' @param ... optional arguments to \code{\link[texreg]{texreg}}.
53
#'
54
#' @details If using \code{custom.coef.map} with more than one model, you should rename
55
#' the parameters in the model objects to ensure that different parameters with the
56
#' same subscript are not conflated by \code{texreg} e.g. \code{beta[1]} could represent age
57
#' in one model and income in another, and \code{texreg} would combine the two if you
58
#' do not rename \code{beta[1]} to more informative names in the model objects.
59
#'
60
#' If \code{mod} is a \code{brmsfit} object or list of \code{brmsfit} objects, note that the
61
#' default \code{brms} names for coefficients are \code{b_Intercept} and \code{b}, so both of
62
#' these should be included in \code{par} if you wish to include the intercept in the
63
#' table.
64
#'
65
#' @return A formatted regression table in LaTeX or HTML format.
66
#'
67
#' @author Rob Williams, \email{jayrobwilliams@gmail.com}
68
#'
69
#' @examples
70
#' \dontshow{.old_wd <- setwd(tempdir())}
71
#' \donttest{
72
#' ## simulating data
73
#' set.seed(123456)
74
#' b0 <- 0.2 # true value for the intercept
75
#' b1 <- 0.5 # true value for first beta
76
#' b2 <- 0.7 # true value for second beta
77
#' n <- 500 # sample size
78
#' X1 <- runif(n, -1, 1)
79
#' X2 <- runif(n, -1, 1)
80
#' Z <- b0 + b1 * X1 + b2 * X2
81
#' pr <- 1 / (1 + exp(-Z)) # inv logit function
82
#' Y <- rbinom(n, 1, pr)
83
#' df <- data.frame(cbind(X1, X2, Y))
84
#' 
85
#' ## formatting the data for jags
86
#' datjags <- as.list(df)
87
#' datjags$N <- length(datjags$Y)
88
#' 
89
#' ## creating jags model
90
#' model <- function()  {
91
#' 
92
#'   for(i in 1:N){
93
#'     Y[i] ~ dbern(p[i])  ## Bernoulli distribution of y_i
94
#'     logit(p[i]) <- mu[i]    ## Logit link function
95
#'     mu[i] <- b[1] +
96
#'       b[2] * X1[i] +
97
#'       b[3] * X2[i]
98
#'   }
99
#' 
100
#'   for(j in 1:3){
101
#'     b[j] ~ dnorm(0, 0.001) ## Use a coefficient vector for simplicity
102
#'   }
103
#' 
104
#' }
105
#' 
106
#' params <- c("b")
107
#' inits1 <- list("b" = rep(0, 3))
108
#' inits2 <- list("b" = rep(0, 3))
109
#' inits <- list(inits1, inits2)
110
#' 
111
#' ## fitting the model with R2jags
112
#' set.seed(123)
113
#' fit <- R2jags::jags(data = datjags, inits = inits,
114
#'                     parameters.to.save = params,
115
#'                     n.chains = 2,
116
#'                     n.iter = 2000, n.burnin = 1000,
117
#'                     model.file = model)
118
#' 
119
#' ## generating regression table with all parameters
120
#' mcmcReg(fit)
121
#' 
122
#' ## generating regression table with only betas and custom coefficent names
123
#' mcmcReg(fit, pars = c('b'), coefnames = c('Variable 1',
124
#'                                           'Variable 2',
125
#'                                           'Variable 3'),
126
#'         regex = TRUE)
127
#' ## generating regression tables with all betas and custom names
128
#' mcmcReg(fit, coefnames = c('Variable 1', 'Variable 2',
129
#'                            'Variable 3', 'deviance'))
130
#' }
131
#' 
132
#' \dontshow{setwd(.old_wd)}
133
#' @export
134
#' 
135
mcmcReg <- function(mod, 
136
                    pars = NULL, 
137
                    pointest = 'mean', 
138
                    ci = .95, 
139
                    hpdi = FALSE,
140
                    sd = FALSE,
141
                    pr = FALSE,
142
                    coefnames = NULL, 
143
                    gof = numeric(0),
144
                    gofnames = character(0),
145
                    format = 'latex', 
146
                    file, 
147
                    regex = FALSE,
148
                    ...) {
149
  
150
  ## pull in unexported functions from other packages
151
  ## other options for future versions might include lifting this and adding authors as copr holders
152 1
  runjags.as.mcmc.list.runjags = getFromNamespace("as.mcmc.list.runjags", "runjags")
153
  
154
  ## if only one model object, coerce to a list
155 1
  if (all(class(mod) != 'list')) mod <- list(mod)
156
  
157
  ## check for heterogeneous model objects
158 1
  if (length(unique(lapply(mod, class))) > 1) stop('More than one object class supplied to argument "mod"')
159
  
160
  ## if only one custom coefficient names vector, coerce to a list
161 1
  if (!is.null(coefnames) & !is.list(coefnames)) coefnames <- list(coefnames)
162
  
163
  ## if only one parameter vector, coerce to a list
164 1
  if (class(pars) != 'list' & !is.null(pars)) pars <- list(pars)
165
  
166
  ## if only one gof statistic scalar or vector, coerce to a list
167 1
  if (class(gof) != 'list') gof <- list(rep(gof, times = length(mod)))
168
  
169
  ## if only one gof statistic name scalar or vector, coerce to a list
170 1
  if (class(gofnames) != 'list') gofnames <- list(gofnames)
171
  
172
  ## extract samples and variable names from jags or rjags objects
173 1
  if (lapply(mod, inherits, what = c('jags', 'rjags'))[[1]]) {
174
    
175
    ## extract posterior samples from list of model objects
176 1
    samps <- lapply(mod, function(x) as.matrix(coda::as.mcmc(x)))
177
    
178
  }
179
  
180
  ## extract samples and variable names from bugs object
181 1
  if (lapply(mod, inherits, what = 'bugs')[[1]]) {
182
    
183
    ## extract posterior samples from list of model objects
184 1
    samps <- lapply(mod, function(x) x$sims.matrix)
185
    
186
  }
187
  
188
  ## extract samples and variable names from runjags object
189 1
  if (lapply(mod, inherits, what = 'runjags')[[1]]) {
190
    
191 1
    samps <- lapply(mod, function(x) as.matrix(runjags.as.mcmc.list.runjags(x)))
192
    
193
  }
194
  
195
  ## extract samples and variable names from remaining objects
196 1
  if (lapply(mod, inherits, what = c("mcmc", "mcmc.list", "stanfit", "stanreg",
197 1
                                     "brmsfit"))[[1]]) {
198
    
199 1
    samps <- lapply(mod, function(x) as.matrix(x))
200
  }
201
  
202
  ## limit samples to supplied parameters
203 1
  if (regex) {
204 1
    samps <- mapply(function(x, y) x[, grepl(x = colnames(x),
205 1
                                            pattern = paste(y, collapse = '|'))],
206 1
                      samps, pars, SIMPLIFY = FALSE)
207 1
  } else if (!is.null(pars)) {
208 1
    samps <- mapply(function(x, y) matrix(x[, y], nrow = nrow(x),
209 1
                                          dimnames = list(NULL, y)),
210 1
                    samps, pars, SIMPLIFY = FALSE)
211
  }
212
  
213
  ## calculate point estimate of posterior density
214 1
  samps_pe <- lapply(samps, function(x) apply(as.matrix(x), 2, get(pointest)))
215
  
216
  ## calculate uncertainty interval for or standard deviation
217 1
  if (sd == TRUE) {
218
    
219 1
    samps_sd <- lapply(samps, function(x) apply(as.matrix(x), 2, sd))
220
    
221 1
  } else if (pr == TRUE) {
222
    
223 1
    samps_sd <- lapply(samps,
224 1
                       function(x) apply(as.matrix(x), 2,
225 1
                                         function(y) mean(sign(y) == sign(mean(y)))))
226
    
227 1
  } else if (hpdi == FALSE) {
228
    
229 1
    samps_ci <- lapply(samps, function(x) apply(as.matrix(x), 2, quantile,
230 1
                                                probs = c(.5 - ci/2, .5 + ci/2)))
231
    
232
  } else {
233
    
234 1
    samps_ci <- lapply(samps, function(x) t(coda::HPDinterval(coda::as.mcmc(x),
235 1
                                                              prob = ci)))
236
    
237
  }
238
  
239
  ## if coefficent names supplied, replace names from model object(s)
240 1
  if (regex & is.null(coefnames)) {
241 1
    coefnames <- mapply(function(x, y) colnames(x)[grepl(x = colnames(x),
242 1
                                                          pattern = paste(y, collapse = '|'))],
243 1
                         samps, pars, SIMPLIFY = FALSE)
244 1
  } else if (is.null(coefnames)) {
245 1
    coefnames <- lapply(samps, colnames)
246
  }
247
  
248
  ##
249 1
  if (length(mod) != length(coefnames)) {
250

251 1
    stop('number of models does not match number of custom coefficient vectors')
252

253
  }
254
  
255
  ## create list of texreg object(s) with point estimates and interval
256 1
  if (sd == TRUE | pr == TRUE) {
257
    
258 1
    tr_list <- mapply(function(v, w, x, y, z) texreg::createTexreg(coef.names = v,
259 1
                                                                   coef = w,
260 1
                                                                   se = x,
261 1
                                                                   gof = y,
262 1
                                                                   gof.names = z),
263 1
                      coefnames, samps_pe, samps_sd, gof, gofnames)
264
    
265
  } else {
266
    
267 1
    tr_list <- mapply(function(v, w, x, y, z) texreg::createTexreg(coef.names = v,
268 1
                                                                   coef = w,
269 1
                                                                   ci.low = x[1, ],
270 1
                                                                   ci.up = x[2, ],
271 1
                                                                   gof = y,
272 1
                                                                   gof.names = z),
273 1
                      coefnames, samps_pe, samps_ci, gof, gofnames)
274
    
275
  }
276
  
277
  ## create LaTeX output
278 1
  if (grepl('tex$', format)) {
279
    
280
    ## create LaTeX code
281 1
    if (sd == TRUE) {
282
      
283 1
      tr <- texreg::texreg(l = tr_list, stars = NULL, ...)
284
      
285 1
    } else if (pr == TRUE) {
286
      
287 1
      tr <- texreg::texreg(l = tr_list, stars = NULL, ...)
288
      
289 1
      tr <- gsub('\\$\\(|\\)\\$', '$', tr)
290
      
291
    } else {
292
      
293 1
      tr <- texreg::texreg(l = tr_list, ...) 
294
      
295
      ## replace confidence w/ credible or highest posterior density in texreg output
296 1
      if (hpdi == FALSE) {
297
        
298 1
        tr <- sub('outside the confidence interval',
299 1
                  paste('outside ', ci * 100 ,'\\\\% credible interval', sep = ''),
300 1
                  tr)
301
        
302
      } else {
303
        
304 1
        tr <- sub('outside the confidence interval',
305 1
                  paste('outside ', ci * 100 ,'\\\\% highest posterior density interval',
306 1
                        sep = ''), tr)
307
        
308
      }
309
    }
310
    
311
    ## return LaTeX code to console or write to file
312 1
    if (missing(file)) {
313
      
314 1
      return(tr)
315
      
316
    } else {
317
      
318
      ## remove newline at start of LaTeX code
319 1
      tr <- sub('^\\n', '', tr)
320
      
321 1
      tex_file <- file(paste(sub('\\.tex$', '', file), 'tex', sep = '.'))
322 1
      writeLines(tr, tex_file, sep = '')
323 1
      close(tex_file)
324
      
325
    }
326
    
327
  }
328
  
329
  ## create HTML output
330 1
  if (format == 'html') {
331
    
332 1
    if (sd == TRUE) {
333
      
334 1
      hr <- texreg::htmlreg(l = tr_list, stars = NULL, ...)
335
      
336 1
    } else if (pr == TRUE) {
337
      
338 1
      hr <- texreg::htmlreg(l = tr_list, stars = NULL, ...)
339
      
340 1
      hr <- gsub('>\\(([0-9]\\.[0-9]{2})\\)<', '>\\1<', hr)
341
      
342
    } else {
343
      
344 1
      hr <- texreg::htmlreg(l = tr_list, ...)
345
      
346
      ## replace confidence w/ credible or highest posterior density in texreg output
347 1
      if (hpdi == FALSE) {
348
        
349 1
        hr <- sub('outside the confidence interval',
350 1
                  paste('outside ', ci * 100, '% credible interval', sep = ''),
351 1
                  hr)
352
        
353
      } else {
354
        
355 1
        hr <- sub('outside the confidence interval',
356 1
                  paste('outside ', ci * 100, '% highest posterior density interval',
357 1
                        sep = ''), hr)
358
        
359
      }
360
      
361
    }
362
    
363
    ## return html code to console or write to file
364 1
    if (missing(file)) {
365
      
366 1
      return(hr)
367
      
368
    } else {
369
      
370 1
      html_file <- file(paste(sub('\\.html$', '', file), 'html', sep = '.'))
371 1
      writeLines(hr, html_file, sep = '')
372 1
      close(html_file)
373
      
374
    }
375
    
376
  }
377
  
378
}

Read our documentation on viewing source code .

Loading