Help with bayesian regression modeling (beginner)

I will take a look to the quoted posts.
Regarding priors, let me see if I understood. As I said, I know that in my field IDIOM forms are usually read slowly than COMP and LEX forms. Also, COMP are read faster than LEX. As this is the “normality” in my field, I should consider a normal prior, let’s say (-1, 1) Question: why you choose -0.5, 0.5 and not -1, 1? Also, TOVIDIOM is generally read slower, thus it is better to put a negative prior like (-1,0)?

On the other side, there are contrastant information about how much fast is the experimental group in comparison to control group, then we cannot be sure about the interaction of TOVGroup. Thus, it is better to not specify a prior for Group and TOVGroup.

###Here the basic model with no priors###

 summary(model)
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: Total_duration_of_whole_fixations ~ Group * TOV + (1 | Participant) + (1 | Media) + (1 | CHR) 
   Data: PHMATPRTLEXGSGC22 (Number of observations: 5236) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Group-Level Effects: 
~CHR (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)    28.19     24.38     1.00    90.24 1.01      407     1023

~Media (Number of levels: 123) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)   129.88      9.84   111.84   150.31 1.00     1031     1628

~Participant (Number of levels: 44) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)   230.74     26.08   183.82   286.16 1.01      428      851

Population-Level Effects: 
                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept            580.93     56.66   469.03   690.55 1.01      274      605
GroupGSPE              4.70     66.24  -120.79   137.98 1.01      250      537
TOVIDIOM              58.97     33.70    -4.13   123.25 1.00      724     1315
TOVLEX                34.58     33.43   -31.08   102.78 1.01      765     1070
GroupGSPE:TOVIDIOM   -18.90     22.28   -61.68    25.91 1.00     3723     2865
GroupGSPE:TOVLEX     -22.89     22.24   -66.14    21.37 1.00     3391     2777

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma   322.42      3.26   316.18   328.86 1.00     4938     2874

Draws were sampled 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).
Messaggio di avvertimento:
There were 4 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup 

#Here the model with normal prior###

However, the code should result as follows:

priors <- get_prior(Total_duration_of_whole_fixations~TOV*Group + (1|Participant)+(1|Media)+(1|CHR), data=PHMATPRTLEXGSGC22, family = gaussian())
> priors
                    prior     class               coef       group resp dpar nlpar lb ub
                   (flat)         b                                                     
                   (flat)         b          GroupGSPE                                  
                   (flat)         b           TOVIDIOM                                  
                   (flat)         b TOVIDIOM:GroupGSPE                                  
                   (flat)         b             TOVLEX                                  
                   (flat)         b   TOVLEX:GroupGSPE                                  
 student_t(3, 513, 317.3) Intercept                                                     
   student_t(3, 0, 317.3)        sd                                                 0   
   student_t(3, 0, 317.3)        sd                            CHR                  0   
   student_t(3, 0, 317.3)        sd          Intercept         CHR                  0   
   student_t(3, 0, 317.3)        sd                          Media                  0   
   student_t(3, 0, 317.3)        sd          Intercept       Media                  0   
   student_t(3, 0, 317.3)        sd                    Participant                  0   
   student_t(3, 0, 317.3)        sd          Intercept Participant                  0   
   student_t(3, 0, 317.3)     sigma                                                 0   
       source
      default
 (vectorized)
 (vectorized)
 (vectorized)
 (vectorized)
 (vectorized)
      default
      default
 (vectorized)
 (vectorized)
 (vectorized)
 (vectorized)
 (vectorized)
 (vectorized)
      default
prior1 <- c(set_prior("normal(-1,1)", class = "b", coef = "TOVIDIOM"))
model <- brm(
  Total_duration_of_whole_fixations ~ Group*TOV + (1 | Participant) + (1 | Media) + (1 | CHR),
  data = PHMATPRTLEXGSGC22,
  family = gaussian(),
  prior = prior1,
  cores = 4,
  chains = 4000
)

As you can see, I didn’t insert sigma because I don’t consider the fixed as slopes. For the same reason, I didn’t insert mean and sd of the dependent variable because the intercept isn’t a slope too (I based my reasoning on the following documentation: 20 Metric Predicted Variable with Multiple Nominal Predictors | Doing Bayesian Data Analysis in brms and the tidyverse (bookdown.org); WAMBS BRMS Tutorial: Popularity Data - Rens van de Schoot.

Then I ran the model:

summary(model)
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: Total_duration_of_whole_fixations ~ Group * TOV + (1 | Participant) + (1 | Media) + (1 | CHR) 
   Data: data1 (Number of observations: 5236) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Group-Level Effects: 
~CHR (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)    27.53     24.89     1.06    90.89 1.01      388     1120

~Media (Number of levels: 123) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)   131.74     10.00   113.72   152.93 1.00      867     1884

~Participant (Number of levels: 44) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)   233.02     26.06   188.43   290.30 1.01      727     1330

Population-Level Effects: 
                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept            609.93     50.64   513.08   709.79 1.01      455     1057
GroupGSPE              0.07     69.67  -135.25   137.94 1.00      411      710
TOVIDIOM              -0.94      0.99    -2.89     1.01 1.00     4992     3108
TOVLEX                 4.68     29.59   -51.45    62.75 1.01      764     1154
GroupGSPE:TOVIDIOM    -7.02     20.79   -48.92    32.16 1.00     2484     2583
GroupGSPE:TOVLEX     -17.41     21.28   -58.89    24.41 1.00     3448     2755

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma   322.47      3.19   316.28   328.57 1.00     5699     2784

Draws were sampled 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).
Messaggio di avvertimento:
There were 5 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup 

I don’t know how to solve the problem of divergence.

> posterior_summary(model)
                                                                      Estimate   Est.Error
b_Intercept                                                       6.099277e+02  50.6354292
b_GroupGSPE                                                       6.684205e-02  69.6677981
b_TOVIDIOM                                                       -9.436583e-01   0.9892072
b_TOVLEX                                                          4.677137e+00  29.5857016
b_GroupGSPE:TOVIDIOM                                             -7.015661e+00  20.7870523
b_GroupGSPE:TOVLEX                                               -1.740856e+01  21.2801373
sd_CHR__Intercept                                                 2.752778e+01  24.8887828
sd_Media__Intercept                                               1.317403e+02  10.0005821
sd_Participant__Intercept                                         2.330224e+02  26.0596789
sigma                                                             3.224679e+02   3.1892878

After that, I thought that TOVIDIOM presents slower RT in comparison to the other levels, but from previous anaylisis I saw that while the differences with TOVCOMP were little, with TOVLEX were wider. Thus, I opted for a (-1,0) prior:

prior1 <- c(set_prior("normal(-1,0)", class = "b", coef = "TOVIDIOM"))

However, the model doesn’t work with this kind of prior:

summary(model)
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: Total_duration_of_whole_fixations ~ Group * TOV + (1 | Participant) + (1 | Media) + (1 | CHR) 
   Data: PHMATPRTLEXGSGC22 (Number of observations: 5236) 

The model does not contain posterior draws.

And even if set put chain=1, another error pops out:

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in 'string', line 57, column 2 to column 37)
Chain 1: 
Chain 1: Initialization between (-2, 2) failed after 100 attempts. 
Chain 1:  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] "Error : Initialization failed."
[1] "error occurred during calling the sampler; sampling not done"

Sorry for the wall of text, I am still trying to understand clearly what I am doing.

Let’s stick with your ###Here the basic model with no priors### first. The get_priors you ran shows you what brms is using as the default priors. So there are no No Priors model in brms, those default priors might make sense given what you know about your data.

I would re-run your basic model but add this

cores=4, chains = 4, iter = 3000,
                  control = list(adapt_delta=0.99, max_treedepth = 12)

inside the model code.

Another helpful thing is to plot your prior distribution in R using. So your normal(-1,0) is likely not to work very well, since it’s centered at -1 with nothing on either side of it. a normal(-0.5,0.5) would be better. That becomes clearer when you plot there. Search around for How to Plot a Normal Distribution in R.

To better see the interactions in your raw data I would plot
library(ggplot2)
ggplot(data1, aes(x=TOV, y=Total_duration_of_whole_fixations, fill=Group, color=Group)) + geom_point(size=1)+geom_jitter(alpha = 0.5, width = 0.3, height = 0.3) +
facet_wrap(~Participant, ncols=4)

Then swap out different things for the facet_wrap.

I am rolling into my usual long weekend here so I will likely be slow to respond until Monday. I’ll write up a bit about how I arrived at my priors and drop that in here.

Hi Ara,
thank you for all the suggestions. Here I put all the computations that I’ve done:

###MODEL###

model <- brm(
  Average_duration_of_whole_fixations ~ Group*TOV + (1 | Participant) + (1 | Media) + (1 | CHR),
  data = dati1
  family = gaussian(),
  cores=4, chains = 4, iter = 3000,
  control = list(adapt_delta=0.99, max_treedepth = 12))

summary(model)
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: Average_duration_of_whole_fixations ~ Group * TOV + (1 | Participant) + (1 | Media) + (1 | CHR) 
   Data: dati1 (Number of observations: 5109) 
  Draws: 4 chains, each with iter = 3000; warmup = 1500; thin = 1;
         total post-warmup draws = 6000

Group-Level Effects: 
~CHR (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     7.28      4.52     0.96    18.70 1.00      923      943

~Media (Number of levels: 123) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     9.09      1.03     7.15    11.19 1.00     2315     3344

~Participant (Number of levels: 44) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)    25.38      2.94    20.28    31.82 1.00     1076     1931

Population-Level Effects: 
                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept            194.49      6.80   180.76   207.52 1.00      776     1782
GroupGSPE              5.13      8.02   -10.46    20.81 1.00      789     1562
TOVIDIOM               3.08      3.07    -2.95     9.08 1.00     3620     4296
TOVLEX                 8.13      3.22     1.78    14.40 1.00     3354     4386
GroupGSPE:TOVIDIOM    -1.30      3.18    -7.64     4.86 1.00     6387     4617
GroupGSPE:TOVLEX      -1.33      3.19    -7.54     5.12 1.00     5781     4978

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    47.14      0.47    46.22    48.08 1.00    10379     4316

Draws were sampled 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).

###FIXED SUMMARY###

summary(model)$fixed
                     Estimate Est.Error   l-95% CI   u-95% CI     Rhat  Bulk_ESS Tail_ESS
Intercept          194.487388  6.803261 180.762917 207.515347 1.004149  775.6492 1782.105
GroupGSPE            5.131045  8.020527 -10.459308  20.813139 1.003339  788.6582 1561.500
TOVIDIOM             3.075234  3.074926  -2.945852   9.079651 1.000753 3620.4500 4296.137
TOVLEX               8.132991  3.220580   1.784159  14.401427 1.000500 3354.4165 4385.592
GroupGSPE:TOVIDIOM  -1.304167  3.184563  -7.638133   4.858713 1.000855 6386.5904 4616.978
GroupGSPE:TOVLEX    -1.332544  3.188886  -7.536324   5.122626 1.000260 5780.9907 4978.009

###RANDOM SUMMARY###

summary(model)$random
$CHR
              Estimate Est.Error  l-95% CI u-95% CI     Rhat Bulk_ESS Tail_ESS
sd(Intercept) 7.283256   4.52259 0.9626881 18.69751 1.001122 922.9301 943.3385

$Media
              Estimate Est.Error l-95% CI u-95% CI     Rhat Bulk_ESS Tail_ESS
sd(Intercept) 9.094831  1.032876 7.147843 11.18662 1.001209 2314.803 3343.962

$Participant
              Estimate Est.Error l-95% CI u-95% CI     Rhat Bulk_ESS Tail_ESS
sd(Intercept) 25.38271  2.938754 20.28169  31.8192 1.001054 1076.188 1931.054

##priors###

prior_summary(model)
                   prior     class               coef       group resp dpar nlpar lb ub
                  (flat)         b                                                     
                  (flat)         b          GroupGSPE                                  
                  (flat)         b GroupGSPE:TOVIDIOM                                  
                  (flat)         b   GroupGSPE:TOVLEX                                  
                  (flat)         b           TOVIDIOM                                  
                  (flat)         b             TOVLEX                                  
 student_t(3, 193, 44.5) Intercept                                                     
   student_t(3, 0, 44.5)        sd                                                 0   
   student_t(3, 0, 44.5)        sd                            CHR                  0   
   student_t(3, 0, 44.5)        sd          Intercept         CHR                  0   
   student_t(3, 0, 44.5)        sd                          Media                  0   
   student_t(3, 0, 44.5)        sd          Intercept       Media                  0   
   student_t(3, 0, 44.5)        sd                    Participant                  0   
   student_t(3, 0, 44.5)        sd          Intercept Participant                  0   
   student_t(3, 0, 44.5)     sigma                                                 0   
       source
      default
 (vectorized)
 (vectorized)
 (vectorized)
 (vectorized)
 (vectorized)
      default
      default
 (vectorized)
 (vectorized)
 (vectorized)
 (vectorized)
 (vectorized)
 (vectorized)

###POSTERIORS###

posterior_summary(model)
                                                                      Estimate   Est.Error
b_Intercept                                                       1.944874e+02  6.80326110
b_GroupGSPE                                                       5.131045e+00  8.02052651
b_TOVIDIOM                                                        3.075234e+00  3.07492567
b_TOVLEX                                                          8.132991e+00  3.22058037
b_GroupGSPE:TOVIDIOM                                             -1.304167e+00  3.18456281
b_GroupGSPE:TOVLEX                                               -1.332544e+00  3.18888605
sd_CHR__Intercept                                                 7.283256e+00  4.52258979
sd_Media__Intercept                                               9.094831e+00  1.03287630
sd_Participant__Intercept                                         2.538271e+01  2.93875359
sigma                                                             4.713869e+01  0.47456831

Just to play a bit, I tried to set other priors:

priors2 <-c(set_prior("normal(-0.5, 0.5)", class="b", coef="TOVIDIOM"), set_prior("normal(-0.5, 0.5)", class="b", coef="GroupGSPE"), set_prior("normal(-0.5, 0.5)", class="b", coef="TOVLEX"))
model <- brm(
+   Average_duration_of_whole_fixations ~ Group*TOV + (1 | Participant) + (1 | Media) +(1 | CHR), data = data1, prior = priors2, family = gaussian(), cores=4, chains = 4, iter = 3000, control = list(adapt_delta=0.99, max_treedepth = 12))
summary(model)
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: Average_duration_of_whole_fixations ~ Group * TOV + (1 | Participant) + (1 | Media) + (1 | CHR) 
   Data: data1 (Number of observations: 5109) 
  Draws: 4 chains, each with iter = 3000; warmup = 1500; thin = 1;
         total post-warmup draws = 6000

Group-Level Effects: 
~CHR (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     6.93      4.00     1.10    16.95 1.00      695      772

~Media (Number of levels: 123) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     9.30      1.04     7.40    11.49 1.00     1802     3037

~Participant (Number of levels: 44) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)    24.95      2.86    20.04    31.30 1.00      809     1538

Population-Level Effects: 
                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept            200.78      5.15   190.15   210.50 1.01      474     1064
GroupGSPE             -0.49      0.49    -1.45     0.51 1.00     5712     4538
TOVIDIOM              -0.51      0.49    -1.46     0.43 1.00     5528     4256
TOVLEX                -0.29      0.49    -1.24     0.67 1.00     6574     4104
GroupGSPE:TOVIDIOM     0.91      2.77    -4.71     6.33 1.00     2987     3622
GroupGSPE:TOVLEX       3.23      2.79    -2.27     8.75 1.00     3057     3736

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    47.14      0.48    46.21    48.10 1.00     6768     4313

Draws were sampled 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).
> posterior_summary(model)
                                                                      Estimate  Est.Error
b_Intercept                                                       2.007770e+02  5.1539194
b_GroupGSPE                                                      -4.867153e-01  0.4947486
b_TOVIDIOM                                                       -5.115936e-01  0.4880417
b_TOVLEX                                                         -2.944787e-01  0.4913180
b_GroupGSPE:TOVIDIOM                                              9.095835e-01  2.7740747
b_GroupGSPE:TOVLEX                                                3.228568e+00  2.7886539
sd_CHR__Intercept                                                 6.933392e+00  4.0030166
sd_Media__Intercept                                               9.300199e+00  1.0378689
sd_Participant__Intercept                                         2.495179e+01  2.8595369
sigma                                                             4.714254e+01  0.4783498

Also, I was wondering. “Group” variables has only two levels. Thus, centering it might make sense. But what about TOV, that present tre levels?

Thank you!

1 Like

Hope your Monday is going well. Looks like things are moving along here.

How many categories are in Group?

The topic of hierarchal models and centered/non-centered gets tangly pretty quickly. @betanalpha has a great post on the topic here
https://betanalpha.github.io/assets/case_studies/hierarchical_modeling.html#1_Modeling_Heterogeneity
that I reference on a regular basis.

As for your priors I would document them and what is your justification for each one. This makes it easier for reviewers or committee members to read and know what you are doing. Typically I would also include all the models I’ve built for a project in a supplemental doc for publication. You will see this called multiverse of models see
https://statmodeling.stat.columbia.edu/2023/02/28/multiverse-r-package/
https://journals.sagepub.com/doi/10.1177/1745691616658637

You may want to dive into the model diagnostics next.

One thing that can help your post be more readable is to add ``` at the start and finished of any code or model summary. I cleaned up your last few posts just to make it a bit easier. Or at the top of the reply you will see </> which will do a similar thing.

Hi Ara,
my Mondy was quite busy, I hope yours was good!
Thank you for the documentation, I will try to have a look tomorrow.

Regarding Group, the levels are two. However, yes maybe we could try some diagnostics for the model. In the meanwhile, I will try to get more experienced about priors.
Thank you also for the scripting suggestions!

1 Like