<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/going-bayes-rstats/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://strengejacke.wordpress.com/2017/08/23/going-bayes-rstats/"> R – Strenge Jacke!</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>Some time ago I started working with Bayesian methods, using the great <a rel="nofollow" href="https://cran.r-project.org/package=rstanarm" target="_blank">rstanarm</a>-package. Beside the fantastic package-vignettes, and books like <em>Statistical Rethinking</em> or <em>Doing Bayesion Data Analysis</em>, I also found the ressources from <a rel="nofollow" href="https://tjmahr.github.io/" target="_blank">Tristan Mahr</a> helpful to both better understand Bayesian analysis and rstanarm. This motivated me to implement tools for Bayesian analysis into my packages, as well.</p>
<p>Due to the latest <a rel="nofollow" href="https://github.com/tidyverse/tidyr" target="_blank">tidyr-update</a>, I had to update some of my packages, in order to make them work again, so – beside some other features – some Bayes-stuff is now avaible in my packages on CRAN.</p>
<h2>Finding shape or location parameters from distributions</h2>
<p>The following functions are included in the <a rel="nofollow" href="https://cran.r-project.org/package=sjstats" target="_blank">sjstats-package</a>. Given some known quantiles or percentiles, or a certain value or ratio and its standard error, the functions <code>find_beta()</code>, <code>find_normal()</code> or <code>find_cauchy()</code> help finding the parameters for a distribution. Taking the example from <a rel="nofollow" href="https://www.johndcook.com/blog/2010/01/31/parameters-from-percentiles/" target="_blank">here</a>, the plot indicates that the mean value for the normal distribution is somewhat above 50. We can find the exact parameters with <code>find_normal()</code>, using the information given in the text:</p>
<pre><code>library(sjstats)
find_normal(x1 = 30, p1 = .1, x2 = 90, p2 = .8)
<span>#> $mean
#> [1] 53.78387
#>
#> $sd
#> [1] 30.48026</span></code></pre>
<h2>High Density Intervals for MCMC samples</h2>
<p>The <code>hdi()</code>-function computes the high density interval for posterior samples. This is nothing special, since there are other packages with such functions as well – however, you can use this function not only on vectors, but also on <code>stanreg</code>-objects (i.e. the results from models fitted with <em>rstanarm</em>). And, if required, you can also transform the HDI-values, e.g. if you need these intervals on an expontiated scale.</p>
<pre><code>library(rstanarm)
fit &lt;- stan_glm(mpg ~ wt + am, data = mtcars, chains = 1)
hdi(fit)
<span>#> term hdi.low hdi.high
#> 1 (Intercept) 32.158505 42.341421
#> 2 wt -6.611984 -4.022419
#> 3 am -2.567573 2.343818
#> 4 sigma 2.564218 3.903652</span>

<span># fit logistic regression model</span>
fit &lt;- stan_glm(
  vs ~ wt + am,
  data = mtcars,
  family = binomial(&quot;logit&quot;),
  chains = 1
)
hdi(fit, prob = .89, trans = exp)
<span>#> term hdi.low hdi.high
#> 1 (Intercept) 4.464230e+02 3.725603e+07
#> 2 wt 6.667981e-03 1.752195e-01
#> 3 am 8.923942e-03 3.747664e-01</span></code></pre>
<h2>Marginal effects for rstanarm-models</h2>
<p>The <a rel="nofollow" href="https://cran.r-project.org/package=ggeffects" target="_blank">ggeffects-package</a> creates tidy data frames of model predictions, which are ready to use with ggplot (though there’s a <code>plot()</code>-method as well). <em>ggeffects</em> supports a wide range of models, and makes it easy to plot marginal effects for specific predictors, includinmg interaction terms. In the past updates, support for more model types was added, for instance <code>polr</code> (pkg MASS), <code>hurdle</code> and <code>zeroinfl</code> (pkg pscl), <code>betareg</code> (pkg betareg), <code>truncreg</code> (pkg truncreg), <code>coxph</code> (pkg survival) and <code>stanreg</code> (pkg rstanarm).</p>
<p><code>ggpredict()</code> is the main function that computes marginal effects. Predictions for <code>stanreg</code>-models are based on the posterior distribution of the linear predictor (<code>posterior_linpred()</code>), mostly for convenience reasons. It is recommended to use the posterior predictive distribution (<code>posterior_predict()</code>) for inference and model checking, and you can do so using the <code>ppd</code>-argument when calling <code>ggpredict()</code>, however, especially for binomial or poisson models, it is harder (and much slower) to compute the „confidence intervals“. That’s why relying on <code>posterior_linpred()</code> is the default for <code>stanreg</code>-models with <code>ggpredict()</code>.</p>
<p>Here is an example with two plots, one without raw data and one including data points:</p>
<pre><code>library(sjmisc)
library(rstanarm)
library(ggeffects)
data(efc)
# make categorical
efc$c161sex &lt;- to_label(efc$c161sex)

# fit model
m &lt;- stan_glm(neg_c_7 ~ c160age + c12hour + c161sex, data = efc)
dat &lt;- ggpredict(m, terms = c(&quot;c12hour&quot;, &quot;c161sex&quot;))
dat
<span>#> # A tibble: 128 x 5
#> x predicted conf.low conf.high group
#> </span></code></pre>