#46 Codecov and expanded tests

Merged Andreas Beger andybega
Missing base report.

Unable to compare commits because the base of the pull request did not upload a coverage report.

Learn more here.


@@ -89,6 +89,8 @@
Loading
89 89
mcmcFDplot <- function(fdfull, 
90 90
                       ROPE = NULL){
91 91
  
92 +
  ROPE <- check_ROPE_argument(ROPE)
93 +
  
92 94
  # convert fdfull to long data frame
93 95
  fd_dat <- tidyr::gather(as.data.frame(fdfull))
94 96
  

@@ -1,155 +1,105 @@
Loading
1 -
#'R function for summarizing MCMC output in a regression-style table
2 -
#'@title Summarize Bayesian MCMC Output
3 -
#'@description R function for summarizing MCMC output in a regression-style table.
4 -
#'@param sims Bayesian model object generated by R2jags, rjags, R2WinBUGS, R2OpenBUGS, 
5 -
#'MCMCpack, rstan, and rstanarm.
6 -
#'@param ci desired level for credible intervals; defaults to c(0.025, 0.975).
7 -
#'@param pars character vector of parameters to be printed; defaults to \code{NULL} 
8 -
#'(all parameters are printed). If not \code{NULL}, the user can either specify the exact names of 
9 -
#'parameters to be printed (e.g. \code{c("alpha", "beta1", "beta2")}) or part of a name 
10 -
#'so that all parameters containing that name will be printed (e.g. \code{"beta"} will print \code{beta1}, \code{beta2}, etc.).
11 -
#'@param Pr print percent of posterior draws with same sign as median; defaults to \code{FALSE}.
12 -
#'@param ROPE defaults to \code{NULL}. If not \code{NULL}, a vector of two values defining the region of 
13 -
#'practical equivalence ("ROPE"); returns \% of posterior draws to the left/right of ROPE. For this quantity 
14 -
#'to be meaningful, all parameters must be on the same scale (e.g. standardized coefficients 
15 -
#'or first differences). See Kruschke (2013, Journal of Experimental 
16 -
#'Psychology 143(2): 573-603) for more on the ROPE.
17 -
#'@references Kruschke, John K. 2013. “Bayesian Estimation Supersedes the T-Test.” Journal of 
18 -
#'Experimental Psychology: General 142 (2): 573–603. https://doi.org/10.1037/a0029146.
19 -
#'@return a data frame containing MCMC summary statistics.
20 -
#'@examples
21 -
#' \donttest{
22 -
#' ## simulating data
23 -
#' set.seed(123456)
24 -
#' b0 <- 0.2 # true value for the intercept
25 -
#' b1 <- 0.5 # true value for first beta
26 -
#' b2 <- 0.7 # true value for second beta
27 -
#' n <- 500 # sample size
28 -
#' X1 <- runif(n, -1, 1)
29 -
#' X2 <- runif(n, -1, 1)
30 -
#' Z <- b0 + b1 * X1 + b2 * X2
31 -
#' pr <- 1 / (1 + exp(-Z)) # inv logit function
32 -
#' Y <- rbinom(n, 1, pr)
33 -
#' data <- data.frame(cbind(X1, X2, Y))
34 -
#' 
35 -
#' ## formatting the data for jags
36 -
#' datjags <- as.list(data)
37 -
#' datjags$N <- length(datjags$Y)
38 -
#' 
39 -
#' ## creating jags model
40 -
#' model <- function()  {
41 -
#' 
42 -
#'   for(i in 1:N){
43 -
#'     Y[i] ~ dbern(p[i])  ## Bernoulli distribution of y_i
44 -
#'     logit(p[i]) <- mu[i]    ## Logit link function
45 -
#'     mu[i] <- b[1] +
46 -
#'       b[2] * X1[i] +
47 -
#'       b[3] * X2[i]
48 -
#'   }
49 -
#' 
50 -
#'   for(j in 1:3){
51 -
#'     b[j] ~ dnorm(0, 0.001) ## Use a coefficient vector for simplicity
52 -
#'   }
53 -
#' 
54 -
#' }
55 -
#' 
56 -
#' params <- c("b")
57 -
#' inits1 <- list("b" = rep(0, 3))
58 -
#' inits2 <- list("b" = rep(0, 3))
59 -
#' inits <- list(inits1, inits2)
60 -
#' 
61 -
#' ## fitting the model with R2jags
62 -
#' set.seed(123)
63 -
#' fit <- R2jags::jags(data = datjags, inits = inits,
64 -
#'                     parameters.to.save = params, n.chains = 2, n.iter = 2000,
65 -
#'                     n.burnin = 1000, model.file = model)
1 +
#' Summarize Bayesian MCMC Output
2 +
#'  
3 +
#' R function for summarizing MCMC output in a regression-style table.
4 +
#'  
5 +
#' @param sims Bayesian model object generated by R2jags, rjags, R2WinBUGS, 
6 +
#'   R2OpenBUGS, MCMCpack, rstan, and rstanarm.
7 +
#' @param ci desired level for credible intervals; defaults to c(0.025, 0.975).
8 +
#' @param pars character vector of parameters to be printed; defaults to \code{NULL} 
9 +
#'   (all parameters are printed). If not \code{NULL}, the user can either specify the exact names of 
10 +
#'   parameters to be printed (e.g. \code{c("alpha", "beta1", "beta2")}) or part of a name 
11 +
#'   so that all parameters containing that name will be printed (e.g. \code{"beta"} will print \code{beta1}, \code{beta2}, etc.).
12 +
#' @param Pr print percent of posterior draws with same sign as median; defaults to \code{FALSE}.
13 +
#' @param ROPE defaults to \code{NULL}. If not \code{NULL}, a vector of two values defining the region of 
14 +
#'   practical equivalence ("ROPE"); returns \% of posterior draws to the left/right of ROPE. For this quantity 
15 +
#'   to be meaningful, all parameters must be on the same scale (e.g. standardized coefficients 
16 +
#'   or first differences). See Kruschke (2013, Journal of Experimental 
17 +
#'   Psychology 143(2): 573-603) for more on the ROPE.
18 +
#'
19 +
#' @references Kruschke, John K. 2013. “Bayesian Estimation Supersedes the T-Test.” Journal of 
20 +
#'   Experimental Psychology: General 142 (2): 573–603. https://doi.org/10.1037/a0029146.
21 +
#'
22 +
#' @return a data frame containing MCMC summary statistics.
23 +
#'
24 +
#' @examples
25 +
#' data("jags_logit")
66 26
#' 
67 27
#' ## printing out table
68 -
#' object <- mcmcTab(fit, 
28 +
#' object <- mcmcTab(jags_logit, 
69 29
#'           ci = c(0.025, 0.975), 
70 30
#'           pars = NULL, 
71 31
#'           Pr = FALSE,
72 32
#'           ROPE = NULL)
73 33
#' object
74 -
#' }
75 -
#'@export
76 -
34 +
#' 
35 +
#' @export
77 36
mcmcTab <- function(sims, 
78 37
                    ci = c(0.025, 0.975), 
79 38
                    pars = NULL, 
80 39
                    Pr = FALSE,
81 40
                    ROPE = NULL) {
82 41
  
83 -
    if(class(sims)[1] == "jags" || class(sims)[1] == "rjags"){
84 -
      sims <- as.matrix(coda::as.mcmc(sims))
85 -
    }
86 -
    if(class(sims)[1] == "bugs"){
87 -
      sims <- sims$sims.matrix
88 -
    }  
89 -
    if(class(sims)[1] == "mcmc"){
90 -
      sims <- as.matrix(sims)
91 -
    }    
92 -
    if(class(sims)[1] == "mcmc.list"){
93 -
      sims <- as.matrix(sims)
94 -
    }      
95 -
    if(class(sims)[1] == "stanreg"){
96 -
      sims <- as.matrix(sims)
97 -
    } 
98 -
    if(class(sims)[1] == "brmsfit"){
99 -
      sims <- as.matrix(sims)
100 -
    } 
101 -
    if(class(sims)[1] == "stanfit"){
102 -
      sims <- as.matrix(sims)
103 -
    }     
104 -
    
105 -
    if(is.null(pars) == TRUE){
106 -
      dat <- sims
107 -
    }
108 -
    
109 -
    if(is.null(pars) == FALSE & length(pars) == 1){
110 -
      dat <- sims[, grepl(x = colnames(sims), pattern = pars)]
111 -
    }
112 -
    
113 -
    if(is.null(pars) == FALSE & length(pars) > 1){
114 -
      dat <- sims[, pars]
115 -
    }
116 -
    
117 -
    dat_wide <- t(dat)
118 -
    
42 +
  if (inherits(sims, what = c("jags", "rjags"))) {
43 +
    sims <- as.matrix(coda::as.mcmc(sims))
44 +
  }
45 +
  if (inherits(sims, what = "bugs")) {
46 +
    sims <- sims$sims.matrix
47 +
  }
48 +
  if (inherits(sims, what = c("mcmc", "mcmc.list", "stanfit", "stanreg",
49 +
                              "brmsfit"))) {
50 +
    sims <- as.matrix(sims)
51 +
  }
52 +
  
53 +
  ROPE <- check_ROPE_argument(ROPE)
54 +
  
55 +
  if(is.null(pars) == TRUE){
56 +
    dat <- sims
57 +
  }
58 +
  
59 +
  if(is.null(pars) == FALSE & length(pars) == 1){
60 +
    dat <- sims[, grepl(x = colnames(sims), pattern = pars)]
61 +
  }
62 +
  
63 +
  if(is.null(pars) == FALSE & length(pars) > 1){
64 +
    dat <- sims[, pars]
65 +
  }
66 +
  
67 +
  dat_wide <- t(dat)
68 +
  
69 +
  mcmctab <- apply(dat_wide, 1, 
70 +
                   function(x) c(Median = round(median(x), digits = 3), # Posterior median
71 +
                                 SD = round(sd(x), digits = 3), # Posterior SD
72 +
                                 Lower = as.numeric(round(quantile(x, probs = ci[1]), digits = 3)), # Lower CI of posterior
73 +
                                 Upper = as.numeric(round(quantile(x, probs = ci[2]), digits = 3)), # Upper CI of posterior
74 +
                                 Pr = round(ifelse(median(x) > 0, length(x[x > 0]) / length(x), length(x[x < 0]) / length(x)), digits = 3) # % of posterior draws with same sign as median
75 +
                   ))
76 +
  
77 +
  if(Pr == FALSE){
119 78
    mcmctab <- apply(dat_wide, 1, 
120 79
                     function(x) c(Median = round(median(x), digits = 3), # Posterior median
121 80
                                   SD = round(sd(x), digits = 3), # Posterior SD
122 81
                                   Lower = as.numeric(round(quantile(x, probs = ci[1]), digits = 3)), # Lower CI of posterior
123 -
                                   Upper = as.numeric(round(quantile(x, probs = ci[2]), digits = 3)), # Upper CI of posterior
124 -
                                   Pr = round(ifelse(median(x) > 0, length(x[x > 0]) / length(x), length(x[x < 0]) / length(x)), digits = 3) # % of posterior draws with same sign as median
125 -
                     ))
126 -
    
127 -
    if(Pr == FALSE){
128 -
      mcmctab <- apply(dat_wide, 1, 
129 -
                       function(x) c(Median = round(median(x), digits = 3), # Posterior median
130 -
                                     SD = round(sd(x), digits = 3), # Posterior SD
131 -
                                     Lower = as.numeric(round(quantile(x, probs = ci[1]), digits = 3)), # Lower CI of posterior
132 -
                                     Upper = as.numeric(round(quantile(x, probs = ci[2]), digits = 3))))
133 -
    }
134 -
    
135 -
    if(!is.null(ROPE)){
136 -
      message("This table contains an estimate for parameter values outside of the region of 
82 +
                                   Upper = as.numeric(round(quantile(x, probs = ci[2]), digits = 3))))
83 +
  }
84 +
  
85 +
  if(!is.null(ROPE)){
86 +
    message("This table contains an estimate for parameter values outside of the region of 
137 87
          practical equivalence (ROPE). For this quantity to be meaningful, all parameters 
138 88
          must be on the same scale (e.g. standardized coefficients or first differences).")
139 -
      
140 -
      mcmctab <- apply(dat_wide, 1, 
141 -
                       function(x) c(Median = round(median(x), digits = 3), # Posterior median
142 -
                                     SD = round(sd(x), digits = 3), # Posterior SD
143 -
                                     Lower = as.numeric(round(quantile(x, probs = ci[1]), digits = 3)), # Lower CI of posterior
144 -
                                     Upper = as.numeric(round(quantile(x, probs = ci[2]), digits = 3)),
145 -
                                     PrOutROPE = round(ifelse(median(x) > 0, length(x[x > ROPE[2]]) / length(x), length(x[x < ROPE[1]]) / length(x)), digits = 3)))
146 -
    }
147 89
    
148 -
    # return(t(mcmctab))
149 -
    out_dat <- data.frame("Variable" = colnames(mcmctab), 
150 -
                          t(mcmctab),
151 -
                          row.names = NULL)
152 -
    
153 -
    return(out_dat)
90 +
    mcmctab <- apply(dat_wide, 1, 
91 +
                     function(x) c(Median = round(median(x), digits = 3), # Posterior median
92 +
                                   SD = round(sd(x), digits = 3), # Posterior SD
93 +
                                   Lower = as.numeric(round(quantile(x, probs = ci[1]), digits = 3)), # Lower CI of posterior
94 +
                                   Upper = as.numeric(round(quantile(x, probs = ci[2]), digits = 3)),
95 +
                                   PrOutROPE = round(ifelse(median(x) > 0, length(x[x > ROPE[2]]) / length(x), length(x[x < ROPE[1]]) / length(x)), digits = 3)))
96 +
  }
97 +
  
98 +
  # return(t(mcmctab))
99 +
  out_dat <- data.frame("Variable" = colnames(mcmctab), 
100 +
                        t(mcmctab),
101 +
                        row.names = NULL)
102 +
  
103 +
  return(out_dat)
154 104
  
155 105
} 

@@ -0,0 +1,12 @@
Loading
1 +
2 +
check_ROPE_argument <- function(x) {
3 +
  
4 +
  if (is.null(x)) return(invisible(x))
5 +
  
6 +
  valid <- is.numeric(x) && length(x)==2 && (x[2] >= x[1])
7 +
  if (!valid) {
8 +
    stop("Invalid ROPE argument; must be a length 2 numeric vector like 'c(0, 1)'")
9 +
  }
10 +
  
11 +
  invisible(x)
12 +
}

Unable to process changes.

No base report to compare against.

Files Coverage
R 81.87%
Project Totals (8 files) 81.87%
Loading