Compute a Bayesian version of R-squared for regression models
Source:R/bayes_R2.R
bayes_R2.brmsfit.RdCompute a Bayesian version of R-squared for regression models
Usage
# S3 method for class 'brmsfit'
bayes_R2(
object,
resp = NULL,
method = NULL,
summary = TRUE,
robust = FALSE,
probs = c(0.025, 0.975),
...
)Arguments
- object
An object of class
brmsfit.- resp
Optional names of response variables. If specified, predictions are performed only for the specified response variables.
- method
Character string specifying how R2 is computed. Either `"model"` to use the model-based variances, or `"residual"` to use the residual-based variances. If `NULL` (default), the model-based R2 is computed where possible, falling back to residual-based R2 for families without current support of model-based variances.
- summary
Should summary statistics be returned instead of the raw values? Default is
TRUE.- robust
If
FALSE(the default) the mean is used as the measure of central tendency and the standard deviation as the measure of variability. IfTRUE, the median and the median absolute deviation (MAD) are applied instead. Only used ifsummaryisTRUE.- probs
The percentiles to be computed by the
quantilefunction. Only used ifsummaryisTRUE.- ...
Further arguments passed to
posterior_epred, which is used in the computation of the R-squared values.
Value
If summary = TRUE, an M x C matrix is returned
(M = number of response variables and c = length(probs) + 2)
containing summary statistics of the Bayesian R-squared values.
If summary = FALSE, the posterior draws of the Bayesian
R-squared values are returned in an S x M matrix (S is the number of draws).
Details
For an introduction to the approach, see Gelman et al. (2019)
and https://github.com/jgabry/bayes_R2/.
For Gaussian and Bernoulli models, bayes_R2 uses model-based residual
variances as proposed in Gelman et al. (2019), with Bernoulli models using
Tjur's pseudo-variance. For Gaussian models with heteroscedastic
sigma, the mean residual variance is used as an approximation (see Tjur
(2009) for discussion on this approximation). For other families,
bayes_R2 warns and falls back to residual-based variances.
References
Andrew Gelman, Ben Goodrich, Jonah Gabry & Aki Vehtari. (2019).
R-squared for Bayesian regression models, The American Statistician,
73(3):307-309. 10.1080/00031305.2018.1549100 (Preprint available at
https://acris.aalto.fi/ws/portalfiles/portal/34206843/bayes_R2_v3.pdf)
Tue Tjur. (2009). Coefficient of determination in logistic regression
models - A new proposal: The coefficient of discrimination, The
American Statistician, 63:366-372. 10.1198/tast.2009.08210
Examples
# \dontrun{
fit <- brm(mpg ~ wt + cyl, data = mtcars)
#> Compiling Stan program...
#> Start sampling
#>
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 6e-06 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 0.019 seconds (Warm-up)
#> Chain 1: 0.017 seconds (Sampling)
#> Chain 1: 0.036 seconds (Total)
#> Chain 1:
#>
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
#> Chain 2:
#> Chain 2: Gradient evaluation took 3e-06 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2:
#> Chain 2:
#> Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 2:
#> Chain 2: Elapsed Time: 0.019 seconds (Warm-up)
#> Chain 2: 0.019 seconds (Sampling)
#> Chain 2: 0.038 seconds (Total)
#> Chain 2:
#>
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
#> Chain 3:
#> Chain 3: Gradient evaluation took 3e-06 seconds
#> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
#> Chain 3: Adjust your expectations accordingly!
#> Chain 3:
#> Chain 3:
#> Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 3:
#> Chain 3: Elapsed Time: 0.021 seconds (Warm-up)
#> Chain 3: 0.021 seconds (Sampling)
#> Chain 3: 0.042 seconds (Total)
#> Chain 3:
#>
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
#> Chain 4:
#> Chain 4: Gradient evaluation took 3e-06 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
#> Chain 4: Adjust your expectations accordingly!
#> Chain 4:
#> Chain 4:
#> Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 4:
#> Chain 4: Elapsed Time: 0.021 seconds (Warm-up)
#> Chain 4: 0.017 seconds (Sampling)
#> Chain 4: 0.038 seconds (Total)
#> Chain 4:
summary(fit)
#> Family: gaussian
#> Links: mu = identity
#> Formula: mpg ~ wt + cyl
#> Data: mtcars (Number of observations: 32)
#> Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#> total post-warmup draws = 4000
#>
#> Regression Coefficients:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept 39.69 1.82 36.07 43.28 1.00 5749 2779
#> wt -3.19 0.79 -4.71 -1.65 1.00 2158 2074
#> cyl -1.51 0.43 -2.36 -0.68 1.00 2248 2015
#>
#> Further Distributional Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 2.69 0.38 2.08 3.56 1.00 2249 2220
#>
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).
bayes_R2(fit)
#> Estimate Est.Error Q2.5 Q97.5
#> R2 0.8047883 0.05285269 0.6779527 0.8837727
# compute R2 with new data
nd <- data.frame(mpg = c(10, 20, 30), wt = c(4, 3, 2), cyl = c(8, 6, 4))
bayes_R2(fit, newdata = nd)
#> Estimate Est.Error Q2.5 Q97.5
#> R2 0.8387443 0.04631315 0.7264823 0.9066032
# }