mcmcRocPrc.R is quite long
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 .