How can I measure the uncertainty interval between two groups of multigroup brms model

Interfaces brms
1

Here is my model

set.seed(1)

n <- 50       # number of reps/group
mu_n <- 0  # means for each group
mu_c1 <- 0.2
mu_c2 <- 0.3
mu_c3 <- 0.3
mu_c4 <- 0.4

d <- tibble(group = rep(c("normal", "covid1", "covid2", "covid3", "covid4"), each = n),
            cases = rep(c (1, 2, 3, 4, 5), each = n)) %>%
   #as.factor(cases) %>%
   mutate(y = c(rnorm(n , mean = mu_n, sd = 1),rnorm(n , mean = mu_c1, sd = 1),rnorm(n , mean = mu_c2, sd = 1),rnorm(n , mean = mu_c3, sd = 1),rnorm(n , mean = mu_c4, sd = 1)))

as.factor(d$cases)

fit <-
brm(data = d,
family = gaussian,
y ~0+ cases,
prior = c(prior(normal(0,2), class = b),
prior(student_t(3,1,1), class = sigma)),
seed =1)

sim_d_and_fit <-
function
(seed, n, probs=c(0.05, 0.95),) {
mu_n <- 0  # means for each group
mu_c1 <- 0.2
mu_c2 <- 0.3
mu_c3 <- 0.3
mu_c4 <- 0.4
set.seed(seed)
d <- tibble(group = rep(c("normal", "covid1", "covid2", "covid3", "covid4"), each = n),
            cases = rep(c (1, 2, 3, 4, 5), each = n)) %>%
   #as.factor(cases) %>%
   mutate(y = c(rnorm(n , mean = mu_n, sd = 1),rnorm(n , mean = mu_c1, sd = 1),rnorm(n , mean = mu_c2, sd = 1),rnorm(n , mean = mu_c3, sd = 1),rnorm(n , mean = mu_c4, sd = 1)))

as.factor(d$cases)

update(fit,
newdata = d,
seed = seed) %>%
fixef(probs = probs) %>%
data.frame() %>%
rownames_to_column("parameter") %>%
filter(parameter =="cases")

n_sim <-100
s1 <-
tibble(seed =1:n_sim) %>%
mutate(b1 = map(seed, sim_d_and_fit, n =60)) %>%
unnest(b1)

My question is how can I measure the UI of the difference between the normal group and covid4 group?? and what does b1 represent here?

2

The coefficients in your brms model will represent the categorical effects via dummy coding. To get to direct statements about the predicted difference between conditions you could use the fitted() function to obtain the expected value for each one, and then generate the contrasts you want from the samples. There are functions to automate this with emmeans. There is also the brmsmargins package to obtain different types of marginal effects, with a detailed tutorial here. The ideal approach will depend on your model and your objective.

1 Like
3

I tired to apply fitted() function

sim_d_and_fit <-
function
(seed, n, probs=c(0.05, 0.95)) {
mu_n <- 0  # means for each group
mu_c1 <- 0.2
mu_c2 <- 0.3
mu_c3 <- 0.3
mu_c4 <- 0.4
set.seed(seed)
d <- tibble(group = rep(c("normal", "covid1", "covid2", "covid3", "covid4"), each = n),
            cases = rep(c (1, 2, 3, 4, 5), each = n)) %>%
   #as.factor(cases) %>%
   mutate(y = c(rnorm(n , mean = mu_n, sd = 1),rnorm(n , mean = mu_c1, sd = 1),rnorm(n , mean = mu_c2, sd = 1),rnorm(n , mean = mu_c3, sd = 1),rnorm(n , mean = mu_c4, sd = 1)))

as.factor(d$cases)

update(fit,
newdata = d,
seed = seed) %>%
fitted()%>%
fixef(probs = probs) %>%
data.frame() %>%
rownames_to_column("parameter") %>%
filter(parameter ==C("normal", "covid4"))# I'm trying to make the model measure the difference between normal and covid 4 group only
}

n_sim <-100
s1 <-
tibble(seed =1:n_sim) %>%
mutate(b1 = map(seed, sim_d_and_fit, n =60)) %>%
unnest(b1)

I had this error

Error in mutate(., b1 = map(seed, sim_d_and_fit, n = 60)) : 
Caused by error in `UseMethod()`:
! no applicable method for 'fixef' applied to an object of class "c('matrix', 'array', 'double', 'numeric')"

Should I do other modifications in the code?

4

@May here the fixef() function is expecting the model object as its first argument but you’re attempting to pass several objects (fixef.brmsfit: Extract Population-Level Estimates in brms: Bayesian Regression Models using 'Stan').

Sorry but I don’t understand the dplyr code so don’t know the exact error here.