-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathhypothesis.brmsfit.Rmd
116 lines (91 loc) · 4.98 KB
/
hypothesis.brmsfit.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
---
title: "``hypothesis.brmsfit`` "
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(brms)
```
Non-Linear Hypothesis Testing
#### Description
Perform non-linear hypothesis testing for all model parameters.
#### Usage
<pre><code>
## S3 method for class 'brmsfit'
hypothesis(x, hypothesis, class = "b", group = "", scope = c("standard", "ranef", "coef"), alpha = 0.05, seed = NULL, ...)
hypothesis(x, ...)
## Default S3 method:
hypothesis(x, hypothesis, alpha = 0.05, ...)
</code></pre>
#### Arguments
* ``x``: An R object. If it is no brmsfit object, it must be coercible to a data.frame.
* ``hypothesis``: A character vector specifying one or more non-linear hypothesis concerning parameters of the model.
* ``class``: A string specifying the class of parameters being tested. Default is "b" for population-level effects. Other typical options are "sd" or "cor". If class = NULL, all parameters can be tested against each other, but have to be specified with their
full name (see also parnames)
* ``group``: Name of a grouping factor to evaluate only group-level effects parameters related to this grouping factor. Ignored if class is not "sd" or "cor".
* ``scope``: Indicates where to look for the variables specified in hypothesis. If "standard", use the full parameter names (subject to the restriction given by class). If "coef" or "ranef" compute the hypothesis for all levels of the grouping factor given in "group", based on the output of coef.brmsfit and ranef.brmsfit, respectively.
* ``alpha``: The alpha-level of the tests (default is 0.05; see 'details' for more information).
* ``seed``: A single numeric value passed to set.seed to make results reproducible.
* ``...``: Currently ignored.
#### Examples
```{r}
## Not run:
## define priors
prior <- c(set_prior("normal(0,2)", class = "b"),
set_prior("student_t(10,0,1)", class = "sigma"),
set_prior("student_t(10,0,1)", class = "sd"))
```
```{r}
## fit a linear mixed effects models
fit <- brm(time ~ age + sex + disease + (1 + age|patient),
data = kidney, family = lognormal(),
prior = prior, sample_prior = TRUE,
control = list(adapt_delta = 0.95))
```
```{r}
## perform two-sided hypothesis testing
(hyp1 <- hypothesis(fit, "sexfemale = age + diseasePKD"))
plot(hyp1)
hypothesis(fit, "exp(age) - 3 = 0", alpha = 0.01)
```
```{r}
## perform one-sided hypothesis testing
hypothesis(fit, "diseasePKD + diseaseGN - 3 < 0")
hypothesis(fit, "age < Intercept",
class = "sd", group = "patient")
```
```{r}
## test the amount of random intercept variance on all variance
h <- paste("sd_patient__Intercept^2 / (sd_patient__Intercept^2 +",
"sd_patient__age^2 + sigma^2) = 0")
(hyp2 <- hypothesis(fit, h, class = NULL))
plot(hyp2)
```
```{r}
## test more than one hypothesis at once
h <- c("diseaseGN = diseaseAN", "2 * diseaseGN - diseasePKD = 0")
(hyp3 <- hypothesis(fit, h))
plot(hyp3, ignore_prior = TRUE)
```
```{r}
## compute hypotheses for all levels of a grouping factor
hypothesis(fit, "age = 0", scope = "coef", group = "patient")
## use the default method
dat <- as.data.frame(fit)
hypothesis(dat, "b_age > 0")
## End(Not run)
```
#### value
A brmshypothesis object.
#### Details
Among others, hypothesis computes an evidence ratio (Evid.Ratio) for each hypothesis. For a directed hypothesis, this is just the posterior probability (Post.Prob) under the hypothesis against its alternative. That is, when the hypothesis if of the form a > b, the evidence ratio is the ratio of the posterior probability of a > b and the posterior probability of a < b. In this example, values greater than one indicate that the evidence in favor of a > b is larger than evidence in favor of a < b.
For an undirected (point) hypothesis, the evidence ratio is a Bayes factor between the hypothesis and its alternative computed via the Savage-Dickey density ratio method. That is the posterior density at the point of interest divided by the prior density at that point. values greater than one indicate that evidence in favor of the point hypothesis has increased after seeing the data.
In order to calculate this Bayes factor, all parameters related to the hypothesis must have proper priors and argument
``sample_prior`` of function brm must be set to TRUE. When interpreting Bayes factors, make sure that your priors are reasonable and carefully chosen, as the result will depend heavily on the priors.
In particular, avoid using default priors. The argument alpha specifies the size of the credible interval (i.e., Bayesian confidence interval).
For instance, if alpha = 0.05 (5%), the credible interval will contain 1 - alpha = 0.95 (95%) of the posterior values. Hence, alpha * 100% of the posterior values will lie outside of the credible interval. Although this allows testing of hypotheses in a similar manner as in the frequentist nullhypothesis testing framework, we strongly argue against using arbitrary cutoffs (e.g., p < .05) to
determine the ’existence’ of an effect.
#### Author(s)
Paul-Christian Buerkner <[email protected]>
See Also
brmshypothesis