sjstats (version 0.12.0)

hdi: Compute high density intervals (HDI) for MCMC samples

Description

hdi() computes the high density interval for values from MCMC samples. rope() calculates the proportion of a posterior distribution that lies within a region of practical equivalence.

Usage

hdi(x, prob = 0.9, trans = NULL, type = c("fixed", "random", "all"))
rope(x, rope, trans = NULL, type = c("fixed", "random", "all"))

Arguments

x

A vector of values from a probability distribution (e.g., posterior probabilities from MCMC sampling), or a stanreg, stanfit, or brmsfit object.

prob

Scalar between 0 and 1, indicating the mass within the credible interval that is to be estimated.

trans

Name of a function or character vector naming a function, used to apply transformations on the returned HDI-values resp. (for rope()) on the values of the posterior distribution, before calculating the rope based on the boundaries given in rope. Note that the values in rope are not transformed.

type

For mixed effects models, specify the type of effects that should be returned. type = "fixed" returns fixed effects only, type = "random" the random effects and type = "all" returns both fixed and random effects.

rope

Vector of length two, indicating the lower and upper limit of a range around zero, which indicates the region of practical equivalence. Values of the posterior distribution within this range are considered as being "practically equivalent to zero".

Value

For hdi(), if x is a vector, returns a vector of length two with the lower and upper limit of the HDI; if x is a stanreg, stanfit or brmsfit object, returns a tibble with lower and upper HDI-limits for each predictor. For rope(), returns the proportion of values from x that are within the boundaries of rope.

Details

Computation for HDI is based on the code from Kruschke 2015, pp. 727f.

References

Kruschke JK. Doing Bayesian Data Analysis: A Tutorial with R, JAGS, and Stan. 2nd edition. Academic Press, 2015

Examples

# NOT RUN {
if (require("rstanarm")) {
  fit <- stan_glm(mpg ~ wt + am, data = mtcars, chains = 1)
  hdi(fit)

  # fit logistic regression model
  fit <- stan_glm(
    vs ~ wt + am,
    data = mtcars,
    family = binomial("logit"),
    chains = 1
  )
  # compute hdi, transform on "odds ratio scale"
  hdi(fit, trans = exp)

  # compute rope, on scale of linear predictor. finds proportion
  # of posterior distribution values between -1 and 1.
  rope(fit, rope = c(-1, 1))

  # compute rope, boundaries as "odds ratios". finds proportion of
  # posterior distribution values, which - after being exponentiated -
  # are between .8 and 1.25 (about -.22 and .22 on linear scale)
  rope(fit, rope = c(.8, 1.25), trans = exp)
}
# }
# NOT RUN {
# }