Skip to content

Commit b8b02db

Browse files
DominiqueMakowskistrengejacke
DominiqueMakowski
and
strengejacke
authored
effective_sample(): returns tail ESS for brms models (#705)
* add * docs * fix * fix * update docs, simplify arguments * do the checks in *insight* * fix test * wordlist * fix test * fix blavaan tests * use posterior for stanmvreg * fix example * update * update rd * rd --------- Co-authored-by: Daniel <[email protected]>
1 parent 624c954 commit b8b02db

34 files changed

+859
-419
lines changed

‎DESCRIPTION

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -104,8 +104,8 @@ Suggests:
104104
parameters,
105105
patchwork,
106106
performance,
107-
quadprog,
108107
posterior,
108+
quadprog,
109109
RcppEigen,
110110
rmarkdown,
111111
rstan,

‎R/describe_posterior.R

Lines changed: 40 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -14,7 +14,7 @@
1414
#' `"rope"`, `"p_map"`, `"p_significance"` (or `"ps"`), `"p_rope"`,
1515
#' `"equivalence_test"` (or `"equitest"`), `"bayesfactor"` (or `"bf"`) or
1616
#' `"all"` to compute all tests. For each "test", the corresponding
17-
#' \pkg{bayestestR} function is called (e.g. [`rope()`] or [`p_direction()`])
17+
#' **bayestestR** function is called (e.g. [`rope()`] or [`p_direction()`])
1818
#' and its results included in the summary output.
1919
#' @param rope_range ROPE's lower and higher bounds. Should be a vector of two
2020
#' values (e.g., `c(-0.1, 0.1)`), `"default"` or a list of numeric vectors of
@@ -49,61 +49,56 @@
4949
#' - [Region of Practical Equivalence (ROPE)](https://easystats.github.io/bayestestR/articles/region_of_practical_equivalence.html)
5050
#' - [Bayes factors](https://easystats.github.io/bayestestR/articles/bayes_factors.html)
5151
#'
52-
#' @examples
52+
#' @examplesIf all(insight::check_if_installed(c("logspline", "rstanarm", "emmeans", "BayesFactor"), quietly = TRUE))
5353
#' library(bayestestR)
5454
#'
55-
#' if (require("logspline")) {
56-
#' x <- rnorm(1000)
57-
#' describe_posterior(x, verbose = FALSE)
58-
#' describe_posterior(x,
59-
#' centrality = "all",
60-
#' dispersion = TRUE,
61-
#' test = "all",
62-
#' verbose = FALSE
63-
#' )
64-
#' describe_posterior(x, ci = c(0.80, 0.90), verbose = FALSE)
55+
#' x <- rnorm(1000)
56+
#' describe_posterior(x, verbose = FALSE)
57+
#' describe_posterior(x,
58+
#' centrality = "all",
59+
#' dispersion = TRUE,
60+
#' test = "all",
61+
#' verbose = FALSE
62+
#' )
63+
#' describe_posterior(x, ci = c(0.80, 0.90), verbose = FALSE)
6564
#'
66-
#' df <- data.frame(replicate(4, rnorm(100)))
67-
#' describe_posterior(df, verbose = FALSE)
68-
#' describe_posterior(
69-
#' df,
70-
#' centrality = "all",
71-
#' dispersion = TRUE,
72-
#' test = "all",
73-
#' verbose = FALSE
74-
#' )
75-
#' describe_posterior(df, ci = c(0.80, 0.90), verbose = FALSE)
65+
#' df <- data.frame(replicate(4, rnorm(100)))
66+
#' describe_posterior(df, verbose = FALSE)
67+
#' describe_posterior(
68+
#' df,
69+
#' centrality = "all",
70+
#' dispersion = TRUE,
71+
#' test = "all",
72+
#' verbose = FALSE
73+
#' )
74+
#' describe_posterior(df, ci = c(0.80, 0.90), verbose = FALSE)
75+
#'
76+
#' df <- data.frame(replicate(4, rnorm(20)))
77+
#' head(reshape_iterations(
78+
#' describe_posterior(df, keep_iterations = TRUE, verbose = FALSE)
79+
#' ))
7680
#'
77-
#' df <- data.frame(replicate(4, rnorm(20)))
78-
#' head(reshape_iterations(
79-
#' describe_posterior(df, keep_iterations = TRUE, verbose = FALSE)
80-
#' ))
81-
#' }
8281
#' \donttest{
8382
#' # rstanarm models
8483
#' # -----------------------------------------------
85-
#' if (require("rstanarm") && require("emmeans")) {
86-
#' model <- suppressWarnings(
87-
#' stan_glm(mpg ~ wt + gear, data = mtcars, chains = 2, iter = 200, refresh = 0)
88-
#' )
89-
#' describe_posterior(model)
90-
#' describe_posterior(model, centrality = "all", dispersion = TRUE, test = "all")
91-
#' describe_posterior(model, ci = c(0.80, 0.90))
92-
#' describe_posterior(model, rope_range = list(c(-10, 5), c(-0.2, 0.2), "default"))
84+
#' model <- suppressWarnings(
85+
#' rstanarm::stan_glm(mpg ~ wt + gear, data = mtcars, chains = 2, iter = 200, refresh = 0)
86+
#' )
87+
#' describe_posterior(model)
88+
#' describe_posterior(model, centrality = "all", dispersion = TRUE, test = "all")
89+
#' describe_posterior(model, ci = c(0.80, 0.90))
90+
#' describe_posterior(model, rope_range = list(c(-10, 5), c(-0.2, 0.2), "default"))
9391
#'
94-
#' # emmeans estimates
95-
#' # -----------------------------------------------
96-
#' describe_posterior(emtrends(model, ~1, "wt"))
97-
#' }
92+
#' # emmeans estimates
93+
#' # -----------------------------------------------
94+
#' describe_posterior(emtrends(model, ~1, "wt"))
9895
#'
9996
#' # BayesFactor objects
10097
#' # -----------------------------------------------
101-
#' if (require("BayesFactor")) {
102-
#' bf <- ttestBF(x = rnorm(100, 1, 1))
103-
#' describe_posterior(bf)
104-
#' describe_posterior(bf, centrality = "all", dispersion = TRUE, test = "all")
105-
#' describe_posterior(bf, ci = c(0.80, 0.90))
106-
#' }
98+
#' bf <- ttestBF(x = rnorm(100, 1, 1))
99+
#' describe_posterior(bf)
100+
#' describe_posterior(bf, centrality = "all", dispersion = TRUE, test = "all")
101+
#' describe_posterior(bf, ci = c(0.80, 0.90))
107102
#' }
108103
#' @export
109104
describe_posterior <- function(posterior, ...) {

‎R/effective_sample.R

Lines changed: 79 additions & 68 deletions
Original file line numberDiff line numberDiff line change
@@ -1,32 +1,67 @@
11
#' Effective Sample Size (ESS)
22
#'
3-
#' This function returns the effective sample size (ESS).
3+
#' Effective Sample Size (ESS) is a measure of how much independent information
4+
#' there is in autocorrelated chains. It is used to assess the quality of MCMC
5+
#' samples. A higher ESS indicates more reliable estimates. For most
6+
#' applications, an effective sample size greater than 1,000 is sufficient for
7+
#' stable estimates (Bürkner, 2017). This function returns the effective sample
8+
#' size (ESS) for various Bayesian model objects. For `brmsfit` objects, the
9+
#' returned ESS corresponds to the bulk-ESS (and the tail-ESS is also returned).
410
#'
511
#' @param model A `stanreg`, `stanfit`, `brmsfit`, `blavaan`, or `MCMCglmm` object.
612
#' @param ... Currently not used.
713
#' @inheritParams hdi
814
#'
15+
#' @inheritSection hdi Model components
16+
#'
917
#' @return A data frame with two columns: Parameter name and effective sample size (ESS).
1018
#'
11-
#' @details **Effective Sample (ESS)** should be as large as possible, altough
19+
#' @details
20+
#' - **Effective Sample (ESS)** should be as large as possible, altough
1221
#' for most applications, an effective sample size greater than 1,000 is
1322
#' sufficient for stable estimates (Bürkner, 2017). The ESS corresponds to the
1423
#' number of independent samples with the same estimation power as the N
1524
#' autocorrelated samples. It is is a measure of \dQuote{how much independent
1625
#' information there is in autocorrelated chains} (*Kruschke 2015, p182-3*).
1726
#'
27+
#' - **Bulk-ESS** is useful as a diagnostic for the sampling efficiency in
28+
#' the bulk of the posterior. It is defined as the effective sample size for
29+
#' rank normalized values using split chains. It can be interpreted as the
30+
#' reliability of indices of central tendency (mean, median, etc.).
31+
#'
32+
#' - **Tail-ESS** is useful as a diagnostic for the sampling efficiency in
33+
#' the tails of the posterior. It is defined as the minimum of the effective
34+
#' sample sizes for 5% and 95% quantiles. It can be interpreted as the
35+
#' reliability of indices that depend on the tails of the distribution (e.g.,
36+
#' credible intervals, tail probabilities, etc.).
37+
#'
1838
#' @references
1939
#' - Kruschke, J. (2014). Doing Bayesian data analysis: A tutorial with R, JAGS,
2040
#' and Stan. Academic Press.
2141
#' - Bürkner, P. C. (2017). brms: An R package for Bayesian multilevel models
2242
#' using Stan. Journal of Statistical Software, 80(1), 1-28
43+
#' - Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., & Bürkner, P.-C.
44+
#' (2021). Rank-normalization, folding, and localization: An improved R-hat
45+
#' for assessing convergence of MCMC. Bayesian Analysis, 16(2), 667-718.
2346
#'
24-
#' @examplesIf require("rstanarm")
47+
#' @examplesIf all(insight::check_if_installed(c("rstanarm", "brms", "posterior"), quietly = TRUE))
2548
#' \donttest{
26-
#' library(rstanarm)
27-
#' model <- suppressWarnings(
28-
#' stan_glm(mpg ~ wt + gear, data = mtcars, chains = 2, iter = 200, refresh = 0)
29-
#' )
49+
#' model <- suppressWarnings(rstanarm::stan_glm(
50+
#' mpg ~ wt + gear,
51+
#' data = mtcars,
52+
#' chains = 2,
53+
#' iter = 200,
54+
#' refresh = 0
55+
#' ))
56+
#' effective_sample(model)
57+
#'
58+
#' model <- suppressWarnings(brms::brm(
59+
#' mpg ~ wt,
60+
#' data = mtcars,
61+
#' chains = 2,
62+
#' iter = 200,
63+
#' refresh = 0
64+
#' ))
3065
#' effective_sample(model)
3166
#' }
3267
#' @export
@@ -50,14 +85,10 @@ effective_sample.default <- function(model, ...) {
5085
#' @rdname effective_sample
5186
#' @export
5287
effective_sample.brmsfit <- function(model,
53-
effects = c("fixed", "random", "all"),
54-
component = c("conditional", "zi", "zero_inflated", "all"),
88+
effects = "fixed",
89+
component = "conditional",
5590
parameters = NULL,
5691
...) {
57-
# check arguments
58-
effects <- match.arg(effects)
59-
component <- match.arg(component)
60-
6192
pars <- insight::find_parameters(
6293
model,
6394
effects = effects,
@@ -66,14 +97,16 @@ effective_sample.brmsfit <- function(model,
6697
flatten = TRUE
6798
)
6899

69-
insight::check_if_installed("rstan")
70-
71-
s <- rstan::summary(model$fit)$summary
72-
s <- subset(s, subset = rownames(s) %in% pars)
100+
insight::check_if_installed("posterior")
101+
idx <- as.data.frame(posterior::summarise_draws(model))
102+
rows_to_keep <- idx$variable %in% pars
103+
# ess_*() functions are defined in:
104+
# https://github.com/stan-dev/posterior/blob/master/R/convergence.R
73105

74106
data.frame(
75-
Parameter = rownames(s),
76-
ESS = round(s[, "n_eff"]),
107+
Parameter = idx$variable[rows_to_keep],
108+
ESS = round(idx[rows_to_keep, "ess_bulk"]),
109+
ESS_tail = round(idx[rows_to_keep, "ess_tail"]),
77110
stringsAsFactors = FALSE,
78111
row.names = NULL
79112
)
@@ -83,57 +116,43 @@ effective_sample.brmsfit <- function(model,
83116
#' @rdname effective_sample
84117
#' @export
85118
effective_sample.stanreg <- function(model,
86-
effects = c("fixed", "random", "all"),
87-
component = c("location", "all", "conditional", "smooth_terms", "sigma", "distributional", "auxiliary"), # nolint
119+
effects = "fixed",
120+
component = "location",
88121
parameters = NULL,
89122
...) {
90-
# check arguments
91-
effects <- match.arg(effects)
92-
component <- match.arg(component)
93-
94-
pars <- insight::find_parameters(
123+
effective_sample.brmsfit(
95124
model,
96125
effects = effects,
97126
component = component,
98127
parameters = parameters,
99-
flatten = TRUE
100-
)
101-
102-
s <- as.data.frame(summary(model))
103-
s <- s[rownames(s) %in% pars, ]
104-
105-
data.frame(
106-
Parameter = rownames(s),
107-
ESS = s[["n_eff"]],
108-
stringsAsFactors = FALSE,
109-
row.names = NULL
128+
...
110129
)
111130
}
112131

113132

114133
#' @export
115134
effective_sample.stanmvreg <- function(model,
116-
effects = c("fixed", "random", "all"),
117-
component = c("location", "all", "conditional", "smooth_terms", "sigma", "distributional", "auxiliary"), # nolint
135+
effects = "fixed",
136+
component = "location",
118137
parameters = NULL,
119138
...) {
120-
# check arguments
121-
effects <- match.arg(effects)
122-
component <- match.arg(component)
123-
124139
pars <- insight::get_parameters(
125140
model,
126141
effects = effects,
127142
component = component,
128143
parameters = parameters
129144
)
130145

131-
s <- as.data.frame(summary(model))
132-
s <- s[rownames(s) %in% colnames(pars), ]
146+
insight::check_if_installed("posterior")
147+
idx <- as.data.frame(posterior::summarise_draws(model))
148+
rows_to_keep <- idx$variable %in% colnames(pars)
149+
# ess_*() functions are defined in:
150+
# https://github.com/stan-dev/posterior/blob/master/R/convergence.R
133151

134152
data.frame(
135-
Parameter = rownames(s),
136-
ESS = s[["n_eff"]],
153+
Parameter = idx$variable[rows_to_keep],
154+
ESS = round(idx[rows_to_keep, "ess_bulk"]),
155+
ESS_tail = round(idx[rows_to_keep, "ess_tail"]),
137156
stringsAsFactors = FALSE,
138157
row.names = NULL
139158
)
@@ -142,18 +161,14 @@ effective_sample.stanmvreg <- function(model,
142161

143162
#' @export
144163
effective_sample.stanfit <- function(model,
145-
effects = c("fixed", "random", "all"),
164+
effects = "fixed",
146165
parameters = NULL,
147166
...) {
148-
# check arguments
149-
effects <- match.arg(effects)
150-
151-
pars <-
152-
insight::get_parameters(
153-
model,
154-
effects = effects,
155-
parameters = parameters
156-
)
167+
pars <- insight::get_parameters(
168+
model,
169+
effects = effects,
170+
parameters = parameters
171+
)
157172

158173
insight::check_if_installed("rstan")
159174

@@ -186,19 +201,15 @@ effective_sample.blavaan <- function(model, parameters = NULL, ...) {
186201

187202
#' @export
188203
effective_sample.MCMCglmm <- function(model,
189-
effects = c("fixed", "random", "all"),
204+
effects = "fixed",
190205
parameters = NULL,
191206
...) {
192-
# check arguments
193-
effects <- match.arg(effects)
194-
195-
pars <-
196-
insight::get_parameters(
197-
model,
198-
effects = effects,
199-
parameters = parameters,
200-
summary = TRUE
201-
)
207+
pars <- insight::get_parameters(
208+
model,
209+
effects = effects,
210+
parameters = parameters,
211+
summary = TRUE
212+
)
202213

203214
s.fixed <- as.data.frame(summary(model)$solutions)
204215
s.random <- as.data.frame(summary(model)$Gcovariances)

0 commit comments

Comments
 (0)