Hi!
I’m a beginner in Bayesian statistics, so I need your help.
I would like to find the main effect using a Bayesian model. The code is the following:
fit <- brm(response ~ 0 + Intercept + contrast + language + contrast:language + (1 | subject), data = test, family = bernoulli("logit"), prior = prior, sample_prior = "yes")
Here, I want to find the evidence regarding the effect of contrast (4 contrasts: heed_hid, head_hid, had_hud, and whod_hood, which is the intercept) on response
What would be the code?
Your help would be much appreciated!
Best,
George
Does summary(fit) give you what you need?
Thank you for your reply.
Main effects are not shown using the summary function (just like in mixed-effects models).
I think it has something to do with hypothesis testing, but I don’t know the exact code. Any ideas?
George
I’m not sure what you mean. For example
out <- brm(
Sepal.Length ~ Sepal.Width + Petal.Length + (1 + Sepal.Width | Species),
data = iris,
cores = 4,
backend = "cmdstanr")
summary(out)
Gives the following output. The population level effects are the main effects. [This output, with its warning about divergent transitions, also indicates that this particular set of posterior samples cannot be reliably used for inference, but this example is just to point you to the main effects.]
> summary(out)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: Sepal.Length ~ Sepal.Width + Petal.Length + (1 + Sepal.Width | Species)
Data: iris (Number of observations: 150)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Group-Level Effects:
~Species (Number of levels: 3)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.86 0.82 0.03 2.99 1.00 926 1182
sd(Sepal.Width) 0.44 0.33 0.06 1.31 1.01 668 447
cor(Intercept,Sepal.Width) 0.04 0.58 -0.93 0.96 1.00 1261 1678
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 1.57 0.59 0.28 2.72 1.00 1568 1107
Sepal.Width 0.42 0.25 -0.10 0.96 1.01 723 592
Petal.Length 0.80 0.06 0.66 0.92 1.00 2625 2146
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.31 0.02 0.28 0.35 1.00 1780 1407
Draws were sampled using sample(hmc). 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).
Warning message:
There were 185 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
So, here, what is the main effect of Petal.Length?
The effect, for example, would be stated as "there is strong evidence for a significant effect of Petal.Length on Sepal.Length
Thus, how can we find the ER for the main effect of an independent variable?
First of all, I realized that the example I whipped up has no interactions. For better clarity, here’s one with both main effects and interactions:
out <- brm(
Sepal.Length ~ Sepal.Width * Petal.Length + (1 + Sepal.Width | Species),
data = iris,
cores = 4,
backend = "cmdstanr")
summary(out)
Yields
> summary(out)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: Sepal.Length ~ Sepal.Width * Petal.Length + (1 + Sepal.Width | Species)
Data: iris (Number of observations: 150)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Group-Level Effects:
~Species (Number of levels: 3)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 1.24 0.98 0.07 3.74 1.01 779 663
sd(Sepal.Width) 0.40 0.34 0.02 1.30 1.01 757 1241
cor(Intercept,Sepal.Width) -0.13 0.57 -0.97 0.93 1.00 1793 2020
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.88 1.19 -1.67 3.02 1.00 841 1414
Sepal.Width 0.64 0.39 -0.16 1.40 1.00 835 1279
Petal.Length 0.99 0.26 0.52 1.52 1.00 627 581
Sepal.Width:Petal.Length -0.06 0.09 -0.24 0.09 1.00 649 787
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.31 0.02 0.27 0.35 1.00 3276 2477
Draws were sampled using sample(hmc). 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).
Warning message:
There were 185 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
The summary provides estimates and 95% credible intervals for each of the three main effects in this model, plus the interaction, all under the heading of “population level effects” as distinct from the group level random effects. I’m not sure from your question what else you want.
Perhaps what @Georgios_Georgiou has in mind is the ‘ANOVA-type’ notion of main effect? This would be a) an omnibus effect of the variable “contrast”, not specific comparisons between the different levels, and b) with the effects of “contrast” estimated across the levels of the other variables that “contrast” is interacting with (e.g., as a mean of the “contrast” effects for each level of “language”).
If that’s what you mean, @Georgios_Georgiou: for b), you can set up the contrasts for “language” to allow that; for a), perhaps a model comparison between models with and without all of the effects of “contrast” (LOO? Bayes factors?).
(although you may not need that at all … the comparisons between different levels of “contrast” are more informative than the omnibus effect of “contrast”)
1 Like
If you want to estimate the main effect instead of the simple effect (treatment coding scheme), switch to a contrast coding scheme. Here’s a link that explains how to do that with brms: 5 The Many Variables & The Spurious Waffles | Statistical rethinking with brms, ggplot2, and the tidyverse: Second edition. If you want to understand the interpretation of different coding schemes, see this book: Chapter 4 Interactions | Learning Statistical Models Through Simulation in R. I don’t know if this is what you mean by main effect but I put the links here for a reference in case it is.