-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathREADME.Rmd
135 lines (91 loc) · 7.11 KB
/
README.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
---
output: github_document
---
<!-- README.md is generated from README.Rmd. Please edit that file -->
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.path = "man/figures/README-",
out.width = "100%"
)
```
# stanova
<!-- badges: start -->
[](https://travis-ci.org/bayesstuff/stanova)
[](https://github.com/bayesstuff/stanova/actions)
[](https://zenodo.org/badge/latestdoi/241858315)
<!-- badges: end -->
The goal of `stanova` is to provide a more relevant and interpretable `summary` for Bayesian models with categorical covariates and possibly interactions and continuous covariates estimated in `Stan`. The core functions are `stanova()` which requires specifying which `rstanarm` function should be called through argument `model_fun` (e.g., `model_fun = glmer` calls `stan_glmer` and allows fitting Bayesian mixed models) and `stanova_brm()` which estimates models using `brms::brm()`.
The issue `stanova` tries to address is that categorical variables with $k$ levels need to be transformed into $k-1$ numerical model coefficients. This poses a problem if a model includes a factor with more than two levels that interacts with another variable. In this case, the most reasonable parameterization of the model is such that the intercept correspond to the (unweighted) grand mean and therefore estimates of the coefficients for the main effects represent average effects (compared to simple effects). In this parameterization, factors with $k$ levels, where $k > 2$ (i.e., more than two levels), cannot be mapped in a 1-to-1 fashion to the $k-1$ coefficients as no such mapping exists. Thus, the estimates of the model coefficients do not represent effects of one factor level, but always pertain to more than one factor level and thus cannot be directly interpreted. In other words, these coefficients should not be looked at. Instead, `stanova` transforms these parameters back such that one gets the information on the factor levels (or combination of factor levels for interactions). The default output shows for each factor level the difference from the intercept which as discussed before corresponds to the (unweighted) grand mean.
Another problem adressed by `stanova` is that for Bayesian models the mapping of factor-levels to model coefficients needs to be done such that the marginal prior for each factor-level is the same. If one uses one of the contrast coding schemes commonly used in frequentist models in which the intercept corresponds to the (unweighted) grand mean the marginal prior differs across factor levels. For example, when using `contr.sum` all but the last factor level are mapped to exactly one model coefficient with positive sign, and the last factor level is mapped with negative sign on all model coefficients. Thus, the marginal prior for the last factor level is more diffuse than for the other factor levels (if $k >2$). `stanova` per default uses the contrast coding scheme suggested by Rouder, Morey, Speckman, and Province (2012) which is such that the marginal prior is the same for all factor levels. When using this contrast, the sum-to-zero constraint that is necessary for the intercept to represent the (unweighted) grand mean is also imposed. This contrast is implemented in the `contr.bayes()` function.
## Installation
For the moment, you can only install the development version from [GitHub](https://github.com/) with:
``` r
# install.packages("remotes")
remotes::install_github("bayesstuff/stanova")
```
At least version `1.5.0` of `emmeans` is needed which can be installed from CRAN:
``` r
install.packages("emmeans")
```
For the `brms` support, in addition to a C++ compiler, at least `brms` version `2.15.0` is needed. This can be installed from CRAN for now:
```r
install.packages("brms")
```
For the `rstanarm` fucntionality, `stanova` attaches the new contrasts to the model instead of to the data if `rstanarm` version `2.21.2`, which is not yet on CRAN, is installed from source:
``` r
Sys.setenv("MAKEFLAGS" = "-j4") ## uses 4 cores during installation
remotes::install_github("stan-dev/rstanarm", build_vignettes = FALSE)
```
## `rstanarm` Example
The most basic example only uses a single factor with three levels to demonstrate the output.
```{r fit1, results='hide', message=FALSE, warning=FALSE}
library(stanova)
data("Machines", package = "MEMSS")
m_machines <- stanova(score ~ Machine + (Machine|Worker),
model_fun = "glmer",
data=Machines, chains = 2,
warmup = 250, iter = 750)
```
The `summary` method of a `stanova` objects first provides some general information about the fitted model (similar to `rstanarm`). It then provides statistics about the intercept. The next block provides estimates for each factor-level of the `Machine` factor. These estimates represent the difference of the factor level against the intercept. For example, for `Machine A`, the estimate is around `-7` suggesting that the mean of this factor level is 7 points below the intercept (i.e., the grand mean).
```{r}
summary(m_machines)
```
If one is not interested in the differences from the factor levels, it is also possible to obtain the estimates marginal means. For this one just needs to set `diff_intercept = FALSE`.
```{r}
summary(m_machines, diff_intercept = FALSE)
```
The key to the output is the `stanova_samples()` function which takes a fitted model objects and returns the posterior samples transformed to represent the difference from the intercept for each factor level (or the marginal means if `diff_intercept = FALSE`). The default output is an `array`, but this can be changed to a `matrix` or `data.frame` with the `return` argument.
```{r}
out_array <- stanova_samples(m_machines)
str(out_array)
```
One can also change which dimension of the `array` represents the chain via the `dimension_chain` argument.
```{r}
out_array2 <- stanova_samples(m_machines, dimension_chain = 2)
str(out_array2)
```
This makes it easy to produce plots via the `bayesplot` package on the level of the differences from the intercept:
```{r, fig.width=8, fig.height=3}
bayesplot::mcmc_trace(out_array2$Machine)
```
## `brms` Example
We can use the same example for `brms`.
```{r fit-brms, results='hide', message=FALSE, warning=FALSE}
library(stanova)
data("Machines", package = "MEMSS")
m_machines_brm <- stanova_brm(score ~ Machine + (Machine|Worker),
data=Machines, chains = 2,
warmup = 250, iter = 750)
```
The `summary` methods works the same. The default shows the difference from the intercept.
```{r}
summary(m_machines_brm)
```
And we can get the estimated marginal means using `diff_intercept = FALSE`.
```{r}
summary(m_machines_brm, diff_intercept = FALSE)
```
# References
Rouder, J. N., Morey, R. D., Speckman, P. L., & Province, J. M. (2012). Default Bayes factors for ANOVA designs. Journal of Mathematical Psychology, 56(5), 356-374. https://doi.org/10.1016/j.jmp.2012.08.001