Predict draws without random effects

Interfaces brms
1

Hi Stan-forum!
I am quite new to Bayesian modeling but have venduered out into multivatiate-hierartichael Bayesian modeling using the brms package in R (version 2.13.5).

I have searched for hours now on my issue but could not find the answer. This might be due to me not knowing the right terms well enough. In that case, I apologize in advance.

I used the brms package to model the effect of three experimental variables on my response variable (replicate devirgence). The variables are selection (yes, no), nutrient (high, low) and time (17 timepoints). I had 3 comparisons within each group which I included in the random-effects (1|comp) in addition to having the variation dependent on time. My formula looks like this:

brm(bf(
replicate divergence ~ time * selection * nutrient + (1 | comp), 
sigma ~time * disturbance)), 
family= gaussian(), 
data=d, iter = 4000, warmup = 2000, prior = prior_normal, control = list(adapt_delta = 0.99))

I now want to plot the model with confidence intervals given only time and selection as visual components. I have managed this using add_fitted_draws() from tidybayes, however this funtion includes all the random effects as well. When I use stat_lineribbon() in the plot the mean line is not smooth like a gaussian model should be.

Does anyone know how to draw fitted_draws while excluding the random effects?

These are my variables/eastimateted parameters:

 get_variables(m_p1_b)
 [1] "b_Intercept"                            "b_sigma_Intercept"                      "b_time"                                 "b_selectionD"                          
 [5] "b_nutrientL"                            "b_time:selectionD"                      "b_time:nutrientL"                       "b_selectionD:nutrientL"                
 [9] "b_time:selectionD:nutrientL"            "b_sigma_time"                           "b_sigma_selectionD"                     "b_sigma_time:selectionD"               
[13] "sd_comparison__Intercept"               "r_comparison[CDH_replica1_2,Intercept]" "r_comparison[CDH_replica1_3,Intercept]" "r_comparison[CDH_replica2_3,Intercept]"
[17] "r_comparison[CDL_replica1_2,Intercept]" "r_comparison[CDL_replica1_3,Intercept]" "r_comparison[CDL_replica2_3,Intercept]" "r_comparison[DCH_replica1_2,Intercept]"
[21] "r_comparison[DCH_replica1_3,Intercept]" "r_comparison[DCH_replica2_3,Intercept]" "r_comparison[DCL_replica1_2,Intercept]" "r_comparison[DCL_replica1_3,Intercept]"
[25] "r_comparison[DCL_replica2_3,Intercept]" "lp__"                                   "accept_stat__"                          "stepsize__"                            
[29] "treedepth__"                            "n_leapfrog__"                           "divergent__"                            "energy__"

And this is the model:

Family: gaussian 
  Links: mu = identity; sigma = log 
Formula: value ~ time * selection * nutrient + (1 | comparison) 
         sigma ~ time * selection
   Data: d (Number of observations: 106) 
Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup samples = 8000

Group-Level Effects: 
~comparison (Number of levels: 12) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.02      0.01     0.00     0.05 1.00     1784     3586

Population-Level Effects: 
                          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept                     0.72      0.03     0.67     0.77 1.00     3422     5429
sigma_Intercept              -2.41      0.10    -2.61    -2.19 1.00     9740     5814
time                         -0.02      0.00    -0.02    -0.01 1.00     4957     4996
selectionD                   -0.24      0.05    -0.34    -0.14 1.00     4943     4995
nutrientL                    -0.04      0.03    -0.11     0.03 1.00     2496     4410
time:selectionD               0.03      0.00     0.02     0.04 1.00     5136     6227
time:nutrientL                0.01      0.00    -0.00     0.01 1.00     5318     4802
selectionD:nutrientL          0.16      0.07     0.03     0.29 1.00     4331     5770
time:selectionD:nutrientL    -0.01      0.01    -0.02     0.00 1.00     5522     5713
sigma_time                    0.06      0.01     0.03     0.09 1.00     4770     5658
sigma_selectionD              0.63      0.15     0.35     0.92 1.00     9900     5355
sigma_time:selectionD        -0.12      0.02    -0.15    -0.08 1.00     4301     5818

Samples were drawn using sampling(NUTS). 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).

Thank you in advance,
Best regards Maddie

2

Hi Maddie, welcome to the Stan forums!

Can you clarify a bit what you mean by this? What do you want on the x-axis of the plot and what do you want on the y-axis? I think if I can clearly understand that then I can hopefully help you figure out the code required to make the plot.

3

Hi Jonah, and thank you for editing my code!

So, currently my code for extracting draws are as follows:

metadata_experiment %>%
  group_by(comparison, nutrient, selection) %>%
  data_grid(time = seq_range(time, n = 9)) %>%
  add_fitted_draws(m_p1_b)

This generates a data frame that I further use ggplot to plot. In my plots I have time on the x axis and the response variable (replicate divergence) on the y axis. In my final plot I want to include both the raw data used to generate the models and the model.

However, as far as I understand, add_fitted_draws() takes into account the group and population effects. But what I am looking for is only the population effects. Do you know a code for that?

4

Yeah I see what you mean. Does adding the argument re_formula = NA to add_fitted_draws() work for you? I think that should tell it to exclude those terms.

1 Like
5

Yes that worked perfectly. Thank you so much.

1 Like