2017 will be the year when social scientists finally decided to diversify their applied statistics toolbox, and stop relying 100% on null hypothesis significance testing (NHST). A very appealing alternative to NHST is Bayesian statistics, which in itself contains many approaches to statistical inference. In this post, I provide an introductory and practical tutorial to Bayesian parameter estimation in the context of comparing two independent groups’ data.

Happy New Year 2017 everybody! 2017 will be the year when social scientists finally decided to diversify their applied statistics toolbox, and stop relying 100% on null hypothesis significance testing (NHST). We now recognize that different scientific questions may require different statistical tools, and are ready to adopt new and innovative methods. A very appealing alternative to NHST is Bayesian statistics, which in itself contains many approaches to statistical inference. In this post, I provide an introductory and practical tutorial to Bayesian parameter estimation in the context of comparing two independent groups’ data.

More specifically, we’ll focus on the t-test. Everyone knows it, everyone uses it. Yet, there are (arguably) better methods for drawing inferences from two independent groups’ metric data (Kruschke 2013):

“When data are interpreted in terms of meaningful parameters in a mathematical description, such as the difference of mean parameters in two groups, it is Bayesian analysis that provides complete information about the credible parameter values. Bayesian analysis is also more intuitive than traditional methods of null hypothesis significance testing (e.g., Dienes, 2011).” (Kruschke 2013)

In that article (*“Bayesian estimation supersedes the t-test”*) Kruschke (2013) provided clear and well-reasoned arguments favoring Bayesian parameter estimation over null hypothesis significance testing in the context of comparing two groups, a situation which is usually dealt with a t-test. It also introduced a “robust” model for comparing two groups, which modeled the data as t-distributed, instead of normal. The article provided R code for running the estimation procedures, which could be downloaded from the author’s website or as an R package.

The R code and programs work well for this specific application (estimating the robust model for one or two groups’ metric data). However, modifying the code to handle more complicated situations is not easy, and the underlying estimation algorithms don’t necessarily scale up to handle more complicated situations. Therefore, in this blog post I’ll introduce easy to use, free, open-source, state-of-the-art computer programs for Bayesian estimation, in the context of comparing two groups’ metric (continuous) data. The programs are available for the R programming language—so make sure you are familiar with R basics (e.g. here). I provide R code for t-tests and Bayesian estimation in R using the R package brms, which provides a concise front-end layer to Stan.

These programs supersede many older Bayesian inference programs because they are easy to use, fast, and are able to handle models with thousands of parameters. Learning to implement basic analyses such as t-tests, and Kruschke’s robust model, with these programs is very useful because you’ll then be able to do Bayesian statistics in practice, and will be prepared to understand and implement more complex models.

Understanding the results of Bayesian estimation requires some knowledge of Bayesian statistics, of course, but since I cannot cover everything in this one post, I refer readers to excellent books on the topic: McElreath (2020), Kruschke (2014), Gelman et al. (2013).

First, I’ll introduce the basic t-test in some detail, and then focus on understanding them as specific instantiations of *linear models*. If that sounds familiar, skip ahead to Bayesian Estimation of the t-test, where I introduce the **brms** package for estimating models using Bayesian methods. Following that, we’ll use “distributional regression” to obtain Bayesian estimates of the unequal variances t-test model. Finally, we’ll learn how to estimate the robust unequal variances model using brms.

We will use the following R packages:

We’ll begin with t-tests, using example data from Kruschke’s paper (p. 577):

“Consider data from two groups of people who take an IQ test. Group 1 (N1=47) consumes a “smart drug,” and Group 2 (N2=42) is a control group that consumes a placebo."

These data are visualized as histograms, below:

These two groups’ IQ scores could be compared with a simple **equal variances t-test** (which you shouldn’t use; Lakens, 2015), also known as Student’s t-test.

```
t.test(IQ ~ Group, data = d, var.equal = T)
```

```
Two Sample t-test
data: IQ by Group
t = -1.5587, df = 87, p-value = 0.1227
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-3.544155 0.428653
sample estimates:
mean in group Control mean in group Treatment
100.3571 101.9149
```

We interpret the t-test in terms of the observed t-value, and whether it exceeds the critical t-value. The critical t-value, in turn, is defined as the extreme \(\alpha / 2\) percentiles of a t-distribution with the given degrees of freedom.

The test results in an observed t-value of 1.56, which is not far enough in the tails of a t-distribution with 87 degrees of freedom to warrant rejecting the null hypothesis (given that we are using \(\alpha\) = .05, which may or may not be an entirely brilliant idea).

Next, we’ll run the more appropriate, **unequal variances t-test** (also known as Welch’s t-test), which R gives by default:

```
t.test(IQ ~ Group, data = d, var.equal = F)
```

```
Welch Two Sample t-test
data: IQ by Group
t = -1.6222, df = 63.039, p-value = 0.1098
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-3.4766863 0.3611848
sample estimates:
mean in group Control mean in group Treatment
100.3571 101.9149
```

Note that while R gives Welch’s t-test by default, SPSS gives both. If you’re using SPSS, make sure to report the Welch’s test results, instead of the equal variances test. Here, the conclusion with respect to rejecting the null hypothesis of equal means is the same. However, notice that the results are numerically different, as they should, because these two t-tests refer to different models.

It is of course up to you, as a researcher, to decide whether you assume equal variances or not. But note that we almost always allow the means to be different (that’s the whole point of the test, really), while many treatments may just as well have an effect on the variances.

The first take-home message from today is that there are actually two t-tests, each associated with a different statistical model. And to make clear what the difference is, we must acquaint ourselves with the models.

We don’t often think of t-*tests* (and ANOVAs) as *models*, but it turns out that they are just linear models disguised as “tests” (see here, here, and here). Recently, there has been a tremendous push for model/parameter estimation, instead of null hypothesis significance testing (Gigerenzer 2004; Cumming 2014; Kruschke 2014), so we will benefit from thinking about t-tests as linear models. Doing so will facilitate seamlessly expanding our models to handle more complicated situations.

The equal variances t-test models metric data with three parameters: Mean for group A, mean for group B, and one shared standard deviation (i.e. the assumption that the standard deviations are equal between the two groups.)

We call the metric outcome variable (IQ scores in our example) \(y_{ik}\), where \(i\) is a subscript indicating the \(i^{th}\) datum, and \(k\) indicates the \(k^{th}\) group. So \(y_{19, 1}\) would be the 19th datum, belonging to group 1. Then we specify that \(y_{ik}\) are normally distributed, \(N(\mu_{ik}, \sigma)\), where \(\mu_{ik}\) indicates the mean of group \(k\), and \(\sigma\) the common standard deviation.

\[y_{ik} \sim N(\mu_{ik}, \sigma^2)\]

Read the formula as “Y is normally distributed with mean \(\mu_{ik}\) (mu), and standard deviation \(\sigma\) (sigma).” Note that the standard deviation \(\sigma\) doesn’t have any subscripts: we assume it is the same for the groups.

The means for groups 0 and 1 are simply \(\mu_0\) and \(\mu_1\), respectively, and their difference (let’s call it \(d\)) is \(d = \mu_0 - \mu_1\). The 95% CI for \(d\) is given in the t-test output, and we can tell that it differs from the one given by Welch’s t-test.

It is unsurprising, then, that if we use a different model (the more appropriate unequal variances model), our inferences may be different. Welch’s t-test is the same as Student’s, except that now we assume (and subsequently estimate) a unique standard deviation \(\sigma_{ik}\) for both groups.

\[y_{ik} \sim N(\mu_{ik}, \sigma_{ik}^2)\]

This model makes a lot of sense, because rarely are we in a situation to *a priori* decide that the variance of scores in Group A is equal to the variance of scores in Group B. If you use the equal variances t-test, you should be prepared to justify and defend this assumption. (Deciding between models—such as between these two t-tests—is one way in which our prior information enters and influences data analysis.)

Armed with this knowledge, we can now see that “conducting a t-test” can be understood as estimating one of these two models. By estimating the model, we obtain t-values, degrees of freedom, and consequently, p-values.

However, for the models described here, it can be easier to think of the t-test as a specific type of the **general linear model**. We can re-write the t-test in an equivalent way, but instead have a specific parameter for the difference in means by writing it as a linear model. The equal variance model can be written as

\[y_{ik} \sim N(\mu_{ik}, \sigma^2)\] \[\mu_{ik} = \beta_0 + \beta_1 Group_{ik}\]

Here, \(\sigma\) is just as before, but we now model the mean with an intercept (control group’s mean, \(\beta_0\)) and the effect of the treatment (\(\beta_1\)). With this model, \(\beta_1\) directly tells us the estimated difference in the two groups. And because it is a parameter in the model, it has an associated standard error, t-value, degrees of freedom, and a p-value. The model can be estimated in R with the following line of code:

```
olsmod <- lm(IQ ~ Group, data = d)
```

The key input here is a model formula, which in R is specified as `outcome ~ predictor`

(`DV ~ IV`

). Using the `lm()`

function, we estimated a linear model predicting `IQ`

from an intercept (automatically included) and a Group parameter. I called this object `olsmod`

for Ordinary Least Squares Model.

R has it’s own model formula syntax, which is well worth learning. The formula in the previous model, `IQ ~ Group`

means that we want to regress IQ on an intercept (which is implicitly included), and group (`Group`

). Besides the formula, we only need to provide the data, which is contained in `d`

.

You can verify that the results are identical to the equal variances t-test above.

```
summary(olsmod)
```

```
Call:
lm(formula = IQ ~ Group, data = d)
Residuals:
Min 1Q Median 3Q Max
-19.9149 -0.9149 0.0851 1.0851 22.0851
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 100.3571 0.7263 138.184 <2e-16 ***
GroupTreatment 1.5578 0.9994 1.559 0.123
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.707 on 87 degrees of freedom
Multiple R-squared: 0.02717, Adjusted R-squared: 0.01599
F-statistic: 2.43 on 1 and 87 DF, p-value: 0.1227
```

Focus on the `GroupTreatment`

row in the estimated coefficients. `Estimate`

is the point estimate (best guess) of the difference in means. `t value`

is the observed t-value (identical to what `t.test()`

reported), and the p-value (`Pr(>|t|)`

) matches as well. The `(Intercept)`

row refers to \(\beta_0\), which is the control group’s mean.

This way of thinking about the model, where we have parameters for one group’s mean, and the effect of the other group, facilitates focusing on the important parameter, the difference, instead of individual means. However, you can of course compute the difference from the means, or the means from one mean and a difference.

Next, I’ll illustrate how to estimate the equal variances t-test using Bayesian methods.

Estimating this model with R, thanks to the Stan and brms teams, is as easy as the linear regression model we ran above. The most important function in the brms package is `brm()`

, for Bayesian Regression Model(ing). The user needs only to input a model formula, just as above, and a data frame that contains the variables specified in the formula. `brm()`

then translates the model into Stan language, and asks Stan to compile the model into C++ and draw samples from the posterior distribution. The result is an R object with the estimated results. We run the model and save the results to `mod_eqvar`

for equal variances model:

```
mod_eqvar <- brm(
IQ ~ Group,
data = d,
cores = 4, # Use 4 cores for parallel processing
file = "iqgroup" # Save results into a file
)
```

The results can be viewed with `summary()`

:

```
summary(mod_eqvar)
```

```
Family: gaussian
Links: mu = identity; sigma = identity
Formula: IQ ~ Group
Data: d (Number of observations: 89)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 100.34 0.72 98.92 101.78 1.00 4208 3116
Group 1.57 0.99 -0.38 3.56 1.00 3781 2917
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 4.71 0.36 4.07 5.47 1.00 3992 2799
Samples were drawn 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).
```

Notice that the model contains three parameters, one of which is the shared standard deviation `sigma`

. Compare the output of the Bayesian model to the one estimated with `lm()`

(OLS):

term | estimate | std.error | brms | Estimate | Est.Error |
---|---|---|---|---|---|

(Intercept) | 100.36 | 0.73 | Intercept | 100.34 | 0.72 |

GroupTreatment | 1.56 | 1.00 | Group | 1.57 | 0.99 |

The point estimates (posterior means in the Bayesian model) and standard errors (SD of the respective posterior distribution) are pretty much identical.

We now know the models behind t-tests, and how to estimate the equal variances t-test using the `t.test()`

, `lm()`

, and `brm()`

functions. We also know how to run Welch’s t-test using `t.test()`

. However, estimating the general linear model version of the unequal variances t-test model is slightly more complicated, because it involves specifying predictors for \(\sigma\), the standard deviation parameter.

We only need a small adjustment to the equal variances model to specify the unequal variances model:

\[y_{ik} \sim N(\mu_{ik}, \sigma_{ik})\] \[\mu_{ik} = \beta_0 + \beta_1 Group_{ik}\]

Notice that we now have subscripts for \(\sigma\), denoting that it varies between groups. In fact, we’ll write out a linear model for the standard deviation parameter.

\[\sigma_{ik} = \gamma_0 + \gamma_1 Group_{ik}\]

The model now includes, instead of a common \(\sigma\), one parameter for Group 0’s standard deviation \(\gamma_0\) (gamma), and one for the effect of Group 1 on the standard deviation \(\gamma_1\), such that group 1’s standard deviation is \(\gamma_0 + \gamma_1\). Therefore, we have 4 free parameters, two means and two standard deviations. (The full specification would include prior distributions for all the parameters, but that topic is outside of the scope of this post.) `brm()`

takes more complicated models by wrapping them inside `bf()`

(short for `brmsformula()`

), which is subsequently entered as the first argument to `brm()`

.

```
uneq_var_frm <- bf(IQ ~ Group, sigma ~ Group)
```

You can see that the formula regresses IQ on Group, such that we’ll have an intercept (implicitly included), and an effect of Group 1. We also model the standard deviation sigma on Group.

```
mod_uneqvar <- brm(
uneq_var_frm,
data = d,
cores = 4,
file = "iqgroup-uv"
)
```

```
Family: gaussian
Links: mu = identity; sigma = log
Formula: IQ ~ Group
sigma ~ Group
Data: d (Number of observations: 89)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 100.35 0.38 99.60 101.05 1.00 5225 2457
sigma_Intercept 0.94 0.12 0.73 1.17 1.00 3689 2586
Group 1.55 0.99 -0.34 3.54 1.00 2343 2608
sigma_Group 0.87 0.16 0.55 1.18 1.00 3699 2299
Samples were drawn 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).
```

The model’s output contains our 4 parameters. `Intercept`

is the mean for group 0, `Group 1`

is the “effect of group 1.” The `sigma_Intercept`

is the standard deviation of Group 0, `sigma_Group`

is the effect of group 1 on the standard deviation (the SD of Group 1 is `sigma_Intercept`

+ `sigma_Group`

). The sigmas are implicitly modeled through a log-link (because they must be positive). To convert them back to the scale of the data, they need to be exponentiated. After taking the exponents of the sigmas, the results look like this:

Parameter | Estimate | Est.Error | Q2.5 | Q97.5 |
---|---|---|---|---|

Intercept | 100.35 | 0.38 | 99.60 | 101.05 |

sigma_Intercept | 2.57 | 0.30 | 2.07 | 3.23 |

Group | 1.55 | 0.99 | -0.34 | 3.54 |

sigma_Group | 2.41 | 0.38 | 1.74 | 3.24 |

Keep in mind that the parameters refer to Group 0’s mean (Intercept) and SD (sigma), and the difference between groups in those values (Group) and (sigma_Group). We now have fully Bayesian estimates of the 4 parameters of the unequal variances t-test model. Finally, let’s move on to the “Robust Bayesian Estimation” model.

Kruschke’s robust model is a comparison of two groups, using five parameters: One mean for each group, one standard deviation for each group, just as in the unequal variances model above. The fifth parameter is a “normality” parameter, \(\nu\) (nu), which means that we are now using a t-distribution to model the data. Using a t-distribution to model the data, instead of a Gaussian, means that the model is less sensitive to extreme values. Here’s what the model looks like:

\[y_{ik} \sim T(\nu, \mu_{ik}, \sigma_{ik})\]

Read the above formula as “Y are random draws from a t-distribution with ‘normality’ parameter \(\nu\), mean \(\mu_{ik}\), and standard deviation \(\sigma_{ik}\).” We have linear models for the means and standard deviations, as above.

This model, as you can see, is almost identical to the unequal variances t-test, but instead uses a t distribution (we assume data are t-distributed), and includes the normality parameter. Using `brm()`

we can still use the unequal variances model, but have to specify the t-distribution. We do this by specifying the `family`

argument to be `student`

(as in Student’s t)

```
Family: student
Links: mu = identity; sigma = log; nu = identity
Formula: IQ ~ Group
sigma ~ Group
Data: d (Number of observations: 89)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 100.52 0.20 100.14 100.91 1.00 5007 2964
sigma_Intercept -0.00 0.19 -0.39 0.37 1.00 3837 3059
Group 1.02 0.43 0.16 1.85 1.00 2380 2415
sigma_Group 0.68 0.25 0.20 1.16 1.00 3549 2746
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
nu 1.84 0.46 1.17 2.95 1.00 2799 1942
Samples were drawn 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).
```

You can compare the results to those in Kruschke’s paper (2013, p.578) to verify that they are nearly identical. There are small discrepancies because of limited number of posterior samples, and because the paper reported posterior modes whereas we focused on means.

Finally, here is how to estimate the model using the original code (Kruschke & Meredith, 2015):

I didn’t actually run that code because after numerous attempts, I was unable to install the rjags package that BEST depends on.

Well, that ended up much longer than what I intended. The aim was both to illustrate the ease of Bayesian modeling in R using brms, and highlight the fact that we can easily move from simple t-tests to more complex (and possibly better) models.

If you’ve followed through, you should be able to conduct Student’s (equal variances) and Welch’s (unequal variances) t-tests in R, and to think about those tests as instantiations of general linear models. Further, you should be able to estimate these models using Bayesian methods.

You should now also be familiar with Kruschke’s robust model for comparing two groups’ metric data, and be able to implement it a few lines of R code. This model found credible differences between two groups, although the frequentist t-tests and models reported p-values well above .05. That should be motivation enough to try robust (Bayesian) models on your own data.

The following software packages were used: R [Version 4.0.3; R Core Team (2020)] and the R-packages *brms* [Version 2.15.0; Bürkner (2017); Bürkner (2018)], *broom* [Version 0.7.5.9000; Robinson, Hayes, and Couch (2021)], *coda* [Version 0.19.4; Plummer et al. (2006)], *dplyr* [Version 1.0.5; Wickham et al. (2021)], *forcats* [Version 0.5.1; Wickham (2021a)], *ggplot2* [Version 3.3.3; Wickham (2016)], *kableExtra* [Version 1.3.4; Zhu (2021)], *knitr* [Version 1.31; Xie (2015)], *purrr* [Version 0.3.4; Henry and Wickham (2020)], *Rcpp* [Version 1.0.6; Eddelbuettel and François (2011); Eddelbuettel and Balamuta (2018)], *readr* [Version 1.4.0; Wickham and Hester (2020)], *scales* [Version 1.1.1; Wickham and Seidel (2020)], *stringr* [Version 1.4.0; Wickham (2019)], *tibble* [Version 3.1.0; Müller and Wickham (2021)], *tidybayes* [Version 2.3.1; Kay (2020)], *tidyr* [Version 1.1.3; Wickham (2021b)], and *tidyverse* [Version 1.3.0; Wickham et al. (2019)].

Bürkner, Paul-Christian. 2017. “brms: An R Package for Bayesian Multilevel Models Using Stan.” *Journal of Statistical Software* 80 (1): 1–28. https://doi.org/10.18637/jss.v080.i01.

———. 2018. “Advanced Bayesian Multilevel Modeling with the R Package brms.” *The R Journal* 10 (1): 395–411. https://doi.org/10.32614/RJ-2018-017.

Cumming, Geoff. 2014. “The New Statistics Why and How.” *Psychological Science* 25 (1): 7–29. https://doi.org/10.1177/0956797613504966.

Eddelbuettel, Dirk, and James Joseph Balamuta. 2018. “Extending extitR with extitC++: A Brief Introduction to extitRcpp.” *The American Statistician* 72 (1): 28–36. https://doi.org/10.1080/00031305.2017.1375990.

Eddelbuettel, Dirk, and Romain François. 2011. “Rcpp: Seamless R and C++ Integration.” *Journal of Statistical Software* 40 (8): 1–18. https://doi.org/10.18637/jss.v040.i08.

Gelman, Andrew, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari, and Donald B. Rubin. 2013. *Bayesian Data Analysis, Third Edition*. Boca Raton: Chapman and Hall/CRC.

Gigerenzer, Gerd. 2004. “Mindless Statistics.” *The Journal of Socio-Economics*, Statistical Significance, 33 (5): 587–606. https://doi.org/10.1016/j.socec.2004.09.033.

Henry, Lionel, and Hadley Wickham. 2020. *Purrr: Functional Programming Tools*. https://CRAN.R-project.org/package=purrr.

Kay, Matthew. 2020. *tidybayes: Tidy Data and Geoms for Bayesian Models*. https://doi.org/10.5281/zenodo.1308151.

Kruschke, John K. 2013. “Bayesian Estimation Supersedes the t Test.” *Journal of Experimental Psychology: General* 142 (2): 573–603. https://doi.org/10.1037/a0029146.

———. 2014. *Doing Bayesian Data Analysis: A Tutorial Introduction with R*. 2nd Edition. Burlington, MA: Academic Press.

McElreath, Richard. 2020. *Statistical Rethinking: A Bayesian Course with Examples in R and Stan*. Second. CRC Texts in Statistical Science. Boca Raton: Taylor and Francis, CRC Press.

Müller, Kirill, and Hadley Wickham. 2021. *Tibble: Simple Data Frames*. https://CRAN.R-project.org/package=tibble.

Plummer, Martyn, Nicky Best, Kate Cowles, and Karen Vines. 2006. “CODA: Convergence Diagnosis and Output Analysis for MCMC.” *R News* 6 (1): 7–11. https://journal.r-project.org/archive/.

R Core Team. 2020. *R: A Language and Environment for Statistical Computing*. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Robinson, David, Alex Hayes, and Simon Couch. 2021. *Broom: Convert Statistical Objects into Tidy Tibbles*.

Wickham, Hadley. 2016. *Ggplot2: Elegant Graphics for Data Analysis*. Springer-Verlag New York. https://ggplot2.tidyverse.org.

———. 2019. *Stringr: Simple, Consistent Wrappers for Common String Operations*. https://CRAN.R-project.org/package=stringr.

———. 2021a. *Forcats: Tools for Working with Categorical Variables (Factors)*. https://CRAN.R-project.org/package=forcats.

———. 2021b. *Tidyr: Tidy Messy Data*. https://CRAN.R-project.org/package=tidyr.

Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” *Journal of Open Source Software* 4 (43): 1686. https://doi.org/10.21105/joss.01686.

Wickham, Hadley, Romain François, Lionel Henry, and Kirill Müller. 2021. *Dplyr: A Grammar of Data Manipulation*. https://CRAN.R-project.org/package=dplyr.

Wickham, Hadley, and Jim Hester. 2020. *Readr: Read Rectangular Text Data*. https://CRAN.R-project.org/package=readr.

Wickham, Hadley, and Dana Seidel. 2020. *Scales: Scale Functions for Visualization*. https://CRAN.R-project.org/package=scales.

Xie, Yihui. 2015. *Dynamic Documents with R and Knitr*. 2nd ed. Boca Raton, Florida: Chapman; Hall/CRC. https://yihui.org/knitr/.

Zhu, Hao. 2021. *kableExtra: Construct Complex Table with ’Kable’ and Pipe Syntax*. https://CRAN.R-project.org/package=kableExtra.

If you see mistakes or want to suggest changes, please create an issue on the source repository.

Text and figures are licensed under Creative Commons Attribution CC BY 4.0. Source code is available at https://github.com/mvuorre/mvuorre.github.io, unless otherwise noted. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".

For attribution, please cite this work as

Vuorre (2017, Jan. 2). Sometimes I R: How to Compare Two Groups with Robust Bayesian Estimation in R. Retrieved from https://mvuorre.github.io/posts/2017-01-02-how-to-compare-two-groups-with-robust-bayesian-estimation-using-r-stan-and-brms/

BibTeX citation

@misc{vuorre2017how, author = {Vuorre, Matti}, title = {Sometimes I R: How to Compare Two Groups with Robust Bayesian Estimation in R}, url = {https://mvuorre.github.io/posts/2017-01-02-how-to-compare-two-groups-with-robust-bayesian-estimation-using-r-stan-and-brms/}, year = {2017} }