<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/anova-freak-and-bayesian-hipster-rstats/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://strengejacke.wordpress.com/2018/03/26/anova-freak-and-bayesian-hipster-rstats/"> R – Strenge Jacke!</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>I’m pleased to announce an update of my <a rel="nofollow" href="https://cran.r-project.org/package=sjstats" target="_blank">sjstats-package</a>. New features are specifically implemented for the Anova and Bayesian statistic and summary functions. Here’s a short overview of what’s new…</p>
<p><span></span></p>
<h1>Anova statistics</h1>
<p>Beside the already implemented functions to calculate <a rel="nofollow" href="https://strengejacke.wordpress.com/2017/07/25/effect-size-statistics-for-anova-tables-rstats/" target="_blank">eta-squared, partial eta-squared and omega-squared</a>, it is now also possible to calculate <em>partial omega-squared</em> statistics for Anova-tables.</p>
<pre><code>library(sjstats)
<span># load sample data</span>
data(efc)

<span># fit linear model</span>
fit &lt;- aov(
  c12hour ~ as.factor(e42dep) + as.factor(c172code) + c160age,
  data = efc
)

omega_sq(fit, partial = TRUE)
<span>#> # A tibble: 3 x 2
#> term partial.omegasq
#> 1 as.factor(e42dep) 0.278 #> 2 as.factor(c172code) 0.00547
#> 3 c160age 0.0649</span></code></pre>
<p>Furthermore, <code>eta_sq()</code> and <code>omega_sq()</code> have a <code>ci.lvl</code>-argument, which – if not <code>NULL</code> – also computes a confidence interval.</p>
<p>For eta-squared, i.e. <code>eta_sq()</code> with <code>partial = FALSE</code>, due to non-symmetry, confidence intervals are based on bootstrap-methods. Confidence intervals for partial omega-squared, i.e. <code>omega_sq()</code> with <code>partial = TRUE</code> – are also based on bootstrapping. In these cases, <code>n</code> indicates the number of bootstrap samples to be drawn to compute the confidence intervals.</p>
<pre><code>eta_sq(fit, partial = TRUE, ci.lvl = .8)
<span>#> # A tibble: 3 x 4
#> term partial.etasq conf.low conf.high
#> 1 as.factor(e42dep) 0.281 0.247 0.310 #> 2 as.factor(c172code) 0.00788 0.00109 0.0160
#> 3 c160age 0.0665 0.0467 0.0885</span></code></pre>
<h1>Bayiasian statistics and summaries</h1>
<p>Some of the summary-functions for Bayesian models were polished in the current update, e.g. <code>hdi()</code>. Let’s fit a (non-sense) model with the great <strong>brms</strong>-package first (you’ll see later why I don’t use <strong>rstanarm</strong> here):</p>
<pre><code>library(lme4)
library(brms)
library(dplyr)
data(sleepstudy)

sleepstudy$mygrp &lt;- sample(1:5, size = 180, replace = TRUE)
sleepstudy %
  group_by(mygrp) %>%
  mutate(mysubgrp = sample(1:30, size = n(), replace = TRUE))

m &lt;- brm(
  Reaction ~ Days + (1 | mygrp / mysubgrp) + (1 | Subject),
  data = sleepstudy
)</code></pre>
<h2>Highest Density Intervals</h2>
<p><code>hdi()</code> calculates by default the 89% HDI for Bayesian regression models:</p>
<pre><code>hdi(m)
<span>#> # A tibble: 3 x 3
#> term hdi.low hdi.high
#> 1 b_Intercept 232.00 268.0
#> 2 b_Days 9.15 11.8
#> 3 sigma 27.20 33.5</span></code></pre>
<p>Now you can calculate the HDI for multiple tresholds (probabilities) at once, using the <code>prob</code>-argument:</p>
<pre><code>hdi(m, prob = c(.5, .8, .95))
<span>#> # A tibble: 3 x 7
#> term hdi.low_0.5 hdi.high_0.5 hdi.low_0.8 hdi.high_0.8 hdi.low_0.95 hdi.high_0.95
#> 1 b_Intercept 244.00 258.0 237.00 265.0 229.00 274.0
#> 2 b_Days 9.97 11.1 9.51 11.6 8.90 12.1
#> 3 sigma 28.90 31.5 27.60 32.5 26.60 33.8</span></code></pre>
<h2>Tidy summary</h2>
<p>A tidy summary, similar to the <code>tidy()</code>-function from the <strong>broom</strong>-package, can be obtained with <code>tidy_stan()</code>. By default, for multilevel models, this function only shows the fixed effects. Use the different options from the <code>type</code>-argument if you want random effects only or both random and fixed effects to be printed.</p>
<pre><code>tidy_stan(m)
<span>#> # A tibble: 3 x 8
#> term estimate std.error hdi.low hdi.high n_eff Rhat mcse
#> 1 b_Intercept 250.0 10.900 233.00 268.0 0.326 1.00 0.312 #> 2 b_Days 10.6 0.835 9.18 11.8 1.000 1.00 0.013
#> 3 sigma 30.2 1.930 27.30 33.4 1.000 1.00 0.202 </span></code></pre>
<h2>Intraclass Correlation Coefficient (ICC)</h2>
<p>The ICC is a useful statistic for random-intercept models, which describes the proportion of variance that is explained by the random effects structure (see <code>?icc</code> for more details). The ICC is based on the variance and correlation components of the model, which you typically get with the <code>VarCorr</code> function. According to the <strong>brms</strong>-package, there’s a good and a bad news: from a developers perspective, the „bad“ news is that <code>VarCorr.brmsfit</code> has recently changed in its internal code structure, so I had to revise my <code>icc()</code>-function to work with <strong>brms</strong>-models again, with success:</p>
<pre><code>icc(m)
<span>#> Bayesian mixed model
#> Family: gaussian (identity)
#> Formula: list() Reaction ~ Days + (1 | mygrp/mysubgrp) + (1 | Subject) list()
#> #> ICC (mygrp): 0.032317
#> ICC (mygrp:mysubgrp): 0.013529
#> ICC (Subject): 0.598439</span></code></pre>
<p>The good news is, that <code>VarCorr.brmsfit</code> has recently changed in its internal code structure. Now it is also possible to get the variance and correlation components <em>for each sample</em> of the model, which allows us to compute the ICC <em>for each sample</em> of the model, which again in turn allows us to compute an ICC including uncertainty intervals very quickly – simply add <code>posterior = TRUE</code> to the <code>icc()</code>-function call:</p>
<pre><code>icc(m, posterior = TRUE)
<span>#> <em># Random Effect Variances and ICC</em>
#> #> Family: gaussian (identity)
#> Formula: Reaction ~ Days + (1 | mygrp/mysubgrp) + (1 | Subject)
#> #> <span>## mygrp</span>
#> ICC: 0.022 (HDI 89%: 0.000-0.104)
#> Between-group: 57.406 (HDI 89%: 0.000-277.588)
#> #> <span>## mygrp:mysubgrp</span>
#> ICC: 0.011 (HDI 89%: 0.000-0.050)
#> Between-group: 28.535 (HDI 89%: 0.000-127.271)
#> #> <span>## Subject</span>
#> ICC: 0.578 (HDI 89%: 0.419-0.737)
#> Between-group: 1449.446 (HDI 89%: 682.362-2413.380)
#> #> <span>## Residuals</span>
#> Within-group: 909.039 (HDI 89%: 722.154-1085.644)</span></code></pre>
<p>By default, the 89% interval around the ICC „point estimate“ is shown in the output, however, the <code>print()</code>-method has an argument to change the interval range to any value, e.g.<br>
<code>print(icc(m, posterior = TRUE), prob = .5, digits = 2)</code>.</p>
<h1>Final words</h1>
<p>The colored output above is not just playing with CSS, the R console output is indeed colored! <img alt="🙂" src="https://i0.wp.com/s0.wp.com/wp-content/mu-plugins/wpcom-smileys/twemoji/2/72x72/1f642.png?w=456&ssl=1" class="wp-smiley" data-recalc-dims="1"></p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://strengejacke.wordpress.com/2018/03/26/anova-freak-and-bayesian-hipster-rstats/"> R – Strenge Jacke!</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/bSnaxvbh1-o" height="1"></a></div>