# Statistics for Bayesian Models

#### 2019-09-08

This vignettes demontrates those functions of the sjstats-package that deal especially with Bayesian models. sjstats provides following functions:

• tidy_stan()
• mediation()

Befor we start, we fit some models, including a mediation-object from the mediation-package, which we use for comparison with brms. The functions work with brmsfit, stanreg and stanfit-objects.

library(sjstats)
library(mediation)
library(brms)

data(jobs)
set.seed(123)

# linear models, for mediation analysis
b1 <- lm(job_seek ~ treat + econ_hard + sex + age, data = jobs)
b2 <- lm(depress2 ~ treat + job_seek + econ_hard + sex + age, data = jobs)

# mediation analysis, for comparison with brms
m1 <- mediate(b1, b2, sims = 1000, treat = "treat", mediator = "job_seek")

m3 <- insight::download_model("brms_zi_2")

## Tidy Summary of Bayesian Models

tidy_stan() is no substitute, but rather a convenient alternative to summary(). The major differences are: tidy_stan()

• focusses on the parameter values (estimates) and gives no information on samples, data, or formula
• calculates the HDI rather than equi-tailed intervals
• separates different model parts, e.g. random from fixed effects, or conditional from zero-inflated models
• and prints everything nicely
tidy_stan(m3)
#>
#> Summary Statistics of Stan-Model
#>
#> # Fixed effects (conditional)
#>
#>    Parameter Estimate Std.Error       HDI(89%)     pd  ESS Rhat MCSE
#>  (Intercept)    -0.84      0.28 [-1.44, -0.29]  96.4%  562 1.01 0.02
#>      persons     0.84      0.09 [ 0.66,  1.06] 100.0%  382 1.01 0.00
#>        child    -1.15      0.09 [-1.29, -0.98] 100.0% 1089 1.00 0.00
#>       camper     0.73      0.09 [ 0.58,  0.89] 100.0% 2724 1.00 0.00
#>
#>
#> # Fixed effects (zero-inflated)
#>
#>    Parameter Estimate Std.Error       HDI(89%)     pd  ESS Rhat MCSE
#>  (Intercept)    -0.64      0.71 [-1.93,  0.52]  83.2%  845 1.00 0.03
#>        child     1.88      0.32 [ 1.40,  2.43] 100.0% 2322 1.00 0.02
#>       camper    -0.83      0.36 [-1.41, -0.24]  99.0% 2277 1.00 0.02
#>
#>
#> # Random effects (conditional) Intercept: persons
#>
#>  Parameter Estimate Std.Error       HDI(89%)    pd ESS Rhat MCSE
#>  persons.1    -0.01      0.10 [-0.38,  0.28] 55.3% 572 1.01 0.00
#>  persons.2     0.02      0.09 [-0.17,  0.30] 61.9% 691 1.01 0.00
#>  persons.3    -0.02      0.08 [-0.26,  0.18] 61.3% 340 1.01 0.00
#>  persons.4     0.00      0.09 [-0.32,  0.33] 51.4% 287 1.01 0.01
#>
#>
#> # Random effects (zero-inflated) Intercept: persons
#>
#>  Parameter Estimate Std.Error       HDI(89%)    pd ESS Rhat MCSE
#>  persons.1     1.28      0.78 [ 0.08,  2.70] 95.7% 811 1.00 0.03
#>  persons.2     0.25      0.68 [-0.90,  1.57] 66.5% 759 1.00 0.03
#>  persons.3    -0.18      0.71 [-1.51,  1.01] 59.7% 871 1.00 0.03
#>  persons.4    -1.29      0.74 [-2.62, -0.01] 94.8% 912 1.00 0.03

Additional statistics in the output are:

• Standard errors (which are actually median absolute deviations, when typical = "median", and standard deviations for typical = "mean").
• The number of effective samples, ESS. Effective Sample Size should be as large as possible, although for most applications, an effective sample size greater than 1,000 is sufficient for stable estimates (Bürkner, 2017). The ESS corresponds to the number of independent samples with the same estimation power as the N autocorrelated samples. It is is a measure of “how much independent information there is in autocorrelated chains” (Kruschke 2015, p182-3).
• Rhat statistics; when Rhat is above 1, it usually indicates that the chain has not yet converged, indicating that the drawn samples might not be trustworthy; drawing more iteration may solve this issue
• Monte Carlo Standard Error (mcse);

By default, the “estimate” is the median of the posterior distribution, but this can be changed with the typical-argument.

tidy_stan(m3, typical = "mean")
#>
#> Summary Statistics of Stan-Model
#>
#> # Fixed effects (conditional)
#>
#>    Parameter Estimate Std.Error       HDI(89%)     pd  ESS Rhat MCSE
#>  (Intercept)    -0.83      0.38 [-1.44, -0.29]  96.4%  562 1.01 0.02
#>      persons     0.84      0.13 [ 0.66,  1.06] 100.0%  382 1.01 0.00
#>        child    -1.15      0.10 [-1.29, -0.98] 100.0% 1089 1.00 0.00
#>       camper     0.73      0.09 [ 0.58,  0.89] 100.0% 2724 1.00 0.00
#>
#>
#> # Fixed effects (zero-inflated)
#>
#>    Parameter Estimate Std.Error       HDI(89%)     pd  ESS Rhat MCSE
#>  (Intercept)    -0.70      0.82 [-1.93,  0.52]  83.2%  845 1.00 0.03
#>        child     1.89      0.32 [ 1.40,  2.43] 100.0% 2322 1.00 0.02
#>       camper    -0.84      0.37 [-1.41, -0.24]  99.0% 2277 1.00 0.02
#>
#>
#> # Random effects (conditional) Intercept: persons
#>
#>  Parameter Estimate Std.Error       HDI(89%)    pd ESS Rhat MCSE
#>  persons.1    -0.04      0.22 [-0.38,  0.28] 55.3% 572 1.01 0.00
#>  persons.2     0.04      0.17 [-0.17,  0.30] 61.9% 691 1.01 0.00
#>  persons.3    -0.03      0.16 [-0.26,  0.18] 61.3% 340 1.01 0.00
#>  persons.4     0.00      0.22 [-0.32,  0.33] 51.4% 287 1.01 0.01
#>
#>
#> # Random effects (zero-inflated) Intercept: persons
#>
#>  Parameter Estimate Std.Error       HDI(89%)    pd ESS Rhat MCSE
#>  persons.1     1.34      0.86 [ 0.08,  2.70] 95.7% 811 1.00 0.03
#>  persons.2     0.31      0.82 [-0.90,  1.57] 66.5% 759 1.00 0.03
#>  persons.3    -0.16      0.83 [-1.51,  1.01] 59.7% 871 1.00 0.03
#>  persons.4    -1.31      0.85 [-2.62, -0.01] 94.8% 912 1.00 0.03

To also show random effects of multilevel models, use the effects-argument.

# printing fixed and random effects of multilevel model
tidy_stan(m3, effects = "all")
#>
#> Summary Statistics of Stan-Model
#>
#> # Fixed effects (conditional)
#>
#>    Parameter Estimate Std.Error       HDI(89%)     pd  ESS Rhat MCSE
#>  (Intercept)    -0.84      0.28 [-1.44, -0.29]  96.4%  562 1.01 0.02
#>      persons     0.84      0.09 [ 0.66,  1.06] 100.0%  382 1.01 0.00
#>        child    -1.15      0.09 [-1.29, -0.98] 100.0% 1089 1.00 0.00
#>       camper     0.73      0.09 [ 0.58,  0.89] 100.0% 2724 1.00 0.00
#>
#>
#> # Fixed effects (zero-inflated)
#>
#>    Parameter Estimate Std.Error       HDI(89%)     pd  ESS Rhat MCSE
#>  (Intercept)    -0.64      0.71 [-1.93,  0.52]  83.2%  845 1.00 0.03
#>        child     1.88      0.32 [ 1.40,  2.43] 100.0% 2322 1.00 0.02
#>       camper    -0.83      0.36 [-1.41, -0.24]  99.0% 2277 1.00 0.02
#>
#>
#> # Random effects (conditional) Intercept: persons
#>
#>  Parameter Estimate Std.Error       HDI(89%)    pd ESS Rhat MCSE
#>  persons.1    -0.01      0.10 [-0.38,  0.28] 55.3% 572 1.01 0.00
#>  persons.2     0.02      0.09 [-0.17,  0.30] 61.9% 691 1.01 0.00
#>  persons.3    -0.02      0.08 [-0.26,  0.18] 61.3% 340 1.01 0.00
#>  persons.4     0.00      0.09 [-0.32,  0.33] 51.4% 287 1.01 0.01
#>
#>
#> # Random effects (zero-inflated) Intercept: persons
#>
#>  Parameter Estimate Std.Error       HDI(89%)    pd ESS Rhat MCSE
#>  persons.1     1.28      0.78 [ 0.08,  2.70] 95.7% 811 1.00 0.03
#>  persons.2     0.25      0.68 [-0.90,  1.57] 66.5% 759 1.00 0.03
#>  persons.3    -0.18      0.71 [-1.51,  1.01] 59.7% 871 1.00 0.03
#>  persons.4    -1.29      0.74 [-2.62, -0.01] 94.8% 912 1.00 0.03

If a model has several components, like zero-inflated models, use the component-argument to decide whether to display all components (default), or only specific parts.

# printing fixed and random effects of multilevel model,
# only zero-inflated model part.
tidy_stan(m3, effects = "all", component = "zero_inflated")
#>
#> Summary Statistics of Stan-Model

By default, 89%-HDI are computed (a convention following McElreath 2015), but other or even multiple HDI can be computed using the prob argument.

# two different HDI for multivariate response model
tidy_stan(m2, prob = c(.5, .95))
#>
#> Summary Statistics of Stan-Model
#>
#> # Fixed effects depress2
#>
#>    Parameter Estimate Std.Error       HDI(50%)       HDI(95%)     pd  ESS Rhat MCSE
#>  (Intercept)     2.21      0.15 [ 2.11,  2.31] [ 1.92,  2.50] 100.0% 6159 1.00 0.00
#>        treat    -0.04      0.04 [-0.06, -0.00] [-0.12,  0.04]  81.7% 7829 1.00 0.00
#>     job_seek    -0.24      0.03 [-0.26, -0.22] [-0.29, -0.18] 100.0% 7601 1.00 0.00
#>    econ_hard     0.15      0.02 [ 0.14,  0.16] [ 0.11,  0.19] 100.0% 8246 1.00 0.00
#>          sex     0.11      0.04 [ 0.08,  0.13] [ 0.02,  0.19]  99.5% 7966 1.00 0.00
#>          age     0.00      0.00 [-0.00,  0.00] [-0.00,  0.00]  61.7% 4317 1.00 0.00
#>
#>
#> # Fixed effects jobseek
#>
#>    Parameter Estimate Std.Error       HDI(50%)       HDI(95%)     pd  ESS Rhat MCSE
#>  (Intercept)     3.67      0.12 [ 3.59,  3.76] [ 3.43,  3.92] 100.0% 6300 1.00 0.00
#>        treat     0.07      0.05 [ 0.03,  0.10] [-0.03,  0.17]  90.3% 8829 1.00 0.00
#>    econ_hard     0.05      0.02 [ 0.03,  0.07] [ 0.01,  0.10]  98.9% 7913 1.00 0.00
#>          sex    -0.01      0.05 [-0.04,  0.03] [-0.10,  0.08]  55.4% 8681 1.00 0.00
#>          age     0.00      0.00 [ 0.00,  0.01] [-0.00,  0.01]  97.5% 5277 1.00 0.00

## Summary of Mediation Analysis

mediation() is another summary function, especially for mediation analysis, i.e. for multivariate response models with casual mediation effects.

Let us recall the models:

f1 <- bf(job_seek ~ treat + econ_hard + sex + age)
f2 <- bf(depress2 ~ treat + job_seek + econ_hard + sex + age)

m2 <- brm(f1 + f2 + set_rescor(FALSE), data = jobs, cores = 4)

Here, treat is the treatment effect, job_seek is the mediator effect, f1 describes the mediator model and f2 describes the outcome model.

mediation() returns a data frame with information on the direct effect (median value of posterior samples from treatment of the outcome model), mediator effect (median value of posterior samples from mediator of the outcome model), indirect effect (median value of the multiplication of the posterior samples from mediator of the outcome model and the posterior samples from treatment of the mediation model) and the total effect (median value of sums of posterior samples used for the direct and indirect effect). The proportion mediated is the indirect effect divided by the total effect.

The simplest call just needs the model-object.

mediation(m2)
#>
#> # Causal Mediation Analysis for Stan Model
#>
#>   Treatment: treat
#>    Mediator: job_seek
#>    Response: depress2
#>
#>                  Estimate    HDI (90%)
#>   Direct effect:    -0.04 [-0.11 0.03]
#> Indirect effect:    -0.02 [-0.04 0.00]
#>    Total effect:    -0.05 [-0.13 0.02]
#>
#> Proportion mediated: 28.14% [-79.57% 135.86%]

Typically, mediation() finds the treatment and mediator variables automatically. If this does not work, use the treatment and mediator arguments to specify the related variable names. For all values, the 90% HDIs are calculated by default. Use prob to calculate a different interval.

Here is a comparison with the mediation package. Note that the summary()-output of the mediation package shows the indirect effect first, followed by the direct effect.

summary(m1)
#>
#> Causal Mediation Analysis
#>
#> Quasi-Bayesian Confidence Intervals
#>
#>                Estimate 95% CI Lower 95% CI Upper p-value
#> ACME            -0.0157      -0.0387         0.01    0.19
#> ADE             -0.0438      -0.1315         0.04    0.35
#> Total Effect    -0.0595      -0.1530         0.02    0.21
#> Prop. Mediated   0.2137      -2.0277         2.70    0.32
#>
#> Sample Size Used: 899
#>
#>
#> Simulations: 1000

mediation(m2, prob = .95)
#>
#> # Causal Mediation Analysis for Stan Model
#>
#>   Treatment: treat
#>    Mediator: job_seek
#>    Response: depress2
#>
#>                  Estimate    HDI (95%)
#>   Direct effect:    -0.04 [-0.12 0.04]
#> Indirect effect:    -0.02 [-0.04 0.01]
#>    Total effect:    -0.05 [-0.15 0.03]
#>
#> Proportion mediated: 28.14% [-178.65% 234.94%]

If you want to calculate mean instead of median values from the posterior samples, use the typical-argument. Furthermore, there is a print()-method, which allows to print more digits.

mediation(m2, typical = "mean", prob = .95) %>% print(digits = 4)
#>
#> # Causal Mediation Analysis for Stan Model
#>
#>   Treatment: treat
#>    Mediator: job_seek
#>    Response: depress2
#>
#>                  Estimate        HDI (95%)
#>   Direct effect:  -0.0395 [-0.1244 0.0450]
#> Indirect effect:  -0.0158 [-0.0400 0.0086]
#>    Total effect:  -0.0553 [-0.1482 0.0302]
#>
#> Proportion mediated: 28.5975% [-178.1953% 235.3902%]

As you can see, the results are similar to what the mediation package produces for non-Bayesian models.

# References

Bürkner, P. C. (2017). brms: An R package for Bayesian multilevel models using Stan. Journal of Statistical Software, 80(1), 1-28

Kruschke, J. (2014). Doing Bayesian data analysis: A tutorial with R, JAGS, and Stan. Academic Press.