hypothesis.brmsfit | R Documentation |
Perform non-linear hypothesis testing for all model parameters.
## S3 method for class 'brmsfit'
hypothesis(
x,
hypothesis,
class = "b",
group = "",
scope = c("standard", "ranef", "coef"),
alpha = 0.05,
robust = FALSE,
seed = NULL,
...
)
hypothesis(x, ...)
## Default S3 method:
hypothesis(x, hypothesis, alpha = 0.05, robust = FALSE, ...)
x |
An |
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 |
group |
Name of a grouping factor to evaluate only group-level effects parameters related to this grouping factor. |
scope |
Indicates where to look for the variables specified in
|
alpha |
The alpha-level of the tests (default is 0.05; see 'Details' for more information). |
robust |
If |
seed |
A single numeric value passed to |
... |
Currently ignored. |
Among others, hypothesis
computes an evidence ratio
(Evid.Ratio
) for each hypothesis. For a one-sided hypothesis, this
is just the posterior probability (Post.Prob
) under the hypothesis
against its alternative. That is, when the hypothesis is 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
two-sided (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 "yes"
.
Otherwise Evid.Ratio
(and Post.Prob
) will be NA
.
Please note that, for technical reasons, we cannot sample from priors of
certain parameters classes. Most notably, these include overall intercept
parameters (prior class "Intercept"
) as well as group-level
coefficients. 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 Evid.Ratio
may sometimes be 0
or Inf
implying very
small or large evidence, respectively, in favor of the tested hypothesis.
For one-sided hypotheses pairs, this basically means that all posterior
draws are on the same side of the value dividing the two hypotheses. In
that sense, instead of 0
or Inf,
you may rather read it as
Evid.Ratio
smaller 1 / S
or greater S
, respectively,
where S
denotes the number of posterior draws used in the
computations.
The argument alpha
specifies the size of the credible interval
(i.e., Bayesian confidence interval). For instance, if we tested a
two-sided hypothesis and set alpha = 0.05
(5%) an, 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 null-hypothesis
testing framework, we strongly argue against using arbitrary cutoffs (e.g.,
p < .05
) to determine the 'existence' of an effect.
A brmshypothesis
object.
Paul-Christian Buerkner paul.buerkner@gmail.com
brmshypothesis
## 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"))
## fit a linear mixed effects models
fit <- brm(time ~ age + sex + disease + (1 + age|patient),
data = kidney, family = lognormal(),
prior = prior, sample_prior = "yes",
control = list(adapt_delta = 0.95))
## perform two-sided hypothesis testing
(hyp1 <- hypothesis(fit, "sexfemale = age + diseasePKD"))
plot(hyp1)
hypothesis(fit, "exp(age) - 3 = 0", alpha = 0.01)
## perform one-sided hypothesis testing
hypothesis(fit, "diseasePKD + diseaseGN - 3 < 0")
hypothesis(fit, "age < Intercept",
class = "sd", group = "patient")
## 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)
## test more than one hypothesis at once
h <- c("diseaseGN = diseaseAN", "2 * diseaseGN - diseasePKD = 0")
(hyp3 <- hypothesis(fit, h))
plot(hyp3, ignore_prior = TRUE)
## 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)
str(dat)
hypothesis(dat, "b_age > 0")
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.