Errors in Modeling in RStudio

Modeling
,
1

I’m trying to do a double hierarchical model in R but coming up with a lot of errors. I’m using the package brms.

The code is shown below. The response variable is Mean.HR, season is a fixed effect, and Date and Turtle are random effects. The sigma at the end was “copied” off of a tutorial that said it would help measure individual residuals.

double_model = bf(Mean.HR ~ season + (1|Date), sigma ~ (1|Turtle))

m1double = brm(double_model,
               data = turtledata.HRmean, 
               warmup = 500, iter = 4000, thin = 2, 
               cores = 2, chains = 2, inits = "random", seed = 12345)

The first error I got was:

Warning in doTryCatch(return(expr), name, parentenv, handler) :
  restarting interrupted promise evaluation

I also got a divergent transition warning, although I have no idea what that means.

There were 187 divergent transitions after warmup.

Lastly, I also got these errosr:

There were 2 chains where the estimated Bayesian Fraction of Missing Information was low

The largest R-hat is 1.06, indicating chains have not mixed.

Since these models are really new to me, I don’t really understand these errors. Please help if all possible, it is greatly appreciated!

2

Welcome to the forum! The first thing I’ll say regarding your errors/warnings is that the Stan team has a lot of really good documentation on these errors and how to potentially resolve them. For example, you might just try typing into Google “rstan AND divergent transitions”. Regardless, I’ll do my best to share my understanding of the things you found after fitting the model (though I’m by no means a technical expert).

I personally have no idea what this warning is, but I do get it from time to time. Generally, I notice that it happens if I try to run a brms model while doing other things on my PC – especially when I have another R session open or when I’m tweaking R/Markdown scripts. I’ve found it to be a harmless warning that something was tried, unable to complete, and then redone.

This is a serious warning that indicates that your model results shouldn’t be trusted. You can learn more about the specifics here, which I encourage because it’s useful to know what Hamiltonian Monte Carlo estimation is doing. Anytime that you see this, you should look at your model and figure out whether there are things you need to tweak. Sometimes, it is sufficient to increase the a control parameter called adapt_delta. In brms, the default for this control is .80, so you can always try increasing this to something like 0.90 or 0.95 to see if the number of divergent transitions is reduced/resolved. In brms you can adjust this value like so:

m1double <- brm(double_model,
                data = turtledata.HRmean,
                warmup = 500, iter = 4000, thin = 2,
                cores = 2, chains = 2, seed = 12345, control = list(adapt_delta = 0.90))

In your case, I suspect that part of the problem is related to inadequate priors that produce a posterior that is difficult to explore. While increasing adapt_delta may help remove some of the divergent transitions, it won’t help with the other problems you have. Adding in better priors may thus help with all of the warnings you got.

This is another warning that your model may not be correctly specified and relates to the efficiency with which the posterior is explored. Like divergent transitions, I recommend reading a little more about the warning and how to probe it here.

With brms, model specification is usually very efficient and rarely is it in need of tweaking unless you are doing some very specific kinds of modeling. As a result, I think that this error might be related to the relatively few number of iterations that you’re keeping. This is where someone with more technical knowledge would be able to give a better rule of thumb, but I generally never do less than 1000 warmup samples and I’ve never run a model that required any thinning. Stan does not have the same degree of autocorrelation in the posterior samples as the standard MCMC programs, so thinning is often not needed (in my experience, and I’d guess in yours since this is a fairly straight forward model).

The link that I provided for the BFMI has a discussion of R-hat directly below that section. The R-hat is one way of looking at whether your chains have mixed, though you should also inspect the trace plots.

In this case, I think this could potentially be related to either the absence of prior and/or the few posterior samples and chains your running. The default number (and the minimum I’ve ever personally run) for Stan is 4 chains. In your particular case, you are running 4500 iterations over 2 chains for a total of 9000 iterations; however, 1000 (500 from each chain) are being used for warmup only, so 3500 iterations are usable posterior samples. My understanding of the thin argument (and I can be wrong here since I’ve never needed to use it) is that you are taking every other of these iterations and so cutting this 3500 post-warmup iterations in half. The result is that your posterior estimates are based on just 1750 posterior samples, which is pretty light.

For comparison, the default settings for brms (and Stan I believe as well) is 2000 iterations over 4 chains. Without any specification (so what would happen if you didn’t specify warmup), half of the iterations are treated as warmup and there is no thinning. This means that each of the 4 chains provide 1000 post-warmup iterations for a total of 4000 posterior samples.

Personal Thoughts
If I were in your situation, I’d try the following trouble-shooting steps in the models. Keep in mind that I have no other information about your data, so these may be totally unhelpful.

  1. Drop the extra arguments and stick to the defaults:
m1double <- brm(double_model,
                data = turtledata.HRmean, seed = 12345)

  1. If you still get errors, then try specifying priors (depending on what the errors are). I believe brms will start with flat priors for regression coefficients and then Student T(3, 0, 10) for intercepts, but I may be wrong about that.

  2. If you have reasonable priors and still get errors, then you can play with the defaults. That depends on the errors, though

1 Like
3

@wgoette wow thank you so much! Those comments really helped me understand this function a lot more. I dropped all the extra arguments but unfortuantely got another error: error occurred during calling the sampler; sampling not done. I looked into the problem and tried some of the suggestions but a lot of it was too complicated for me to understand (this is my very first time with the stan programs).

Any ideas on how to fix that?

4

I’m not sure off hand whether I’ve ever run into that issue before. Just a quick Google of the issue makes it seem like this is a more technical error from rstan. Any chance you had a system update, updated R, or updated any Stan-related packages since you last ran the model? If so, you may try re-installing all the rstan stuff to see if that fixes any issues.

Did the model compile really quickly when you tried to run the model? I think maybe sometimes rstan will get tripped up if there is already a compiled model, so you should always remove the old fits before running the model again (or using update(...) instead of a fresh brm(...) call).

If those things don’t work, then it may be useful to provide your session information (e.g., what versions of packages you have, what R version you’re running, whether or not you’re running in RStudio, your computer’s OS, etc) so that someone else with more knowledge of why the error is occurring can help troubleshoot.

5

Since I last ran it, I haven’t done any updates. I’m actually currently updating my computer as some places say that might help.

And the model takes about the same time to run as normal. It lists all the output and has that error attached. I do summary() of the model and it tells me:

Family: gaussian 
  Links: mu = identity; sigma = log 
Formula: Mean.HR ~ season + (1 | Date) 
         sigma ~ (1 | Turtle)
   Data: turtledata.HRmean (Number of observations: 119) 

The model does not contain posterior samples.

I’m not sure what you mean by removing the old fits. Instead of using brm() do I use update()? I tried that and it didn’t work.

But, for further reference, here is some information:

R version 3.6.2
I’m in RStudio: RStudio Desktop 1.4.1717
Mac: macOS Big Sur: 11.2.3 (currently updating to 11.2.4)
brms version: 2.15.0

6

Yeah, I might not be the one to know what’s going on here. I do seem to recall some threads on this forum about the Mac OS setup for rstan – I’ve never read them in-depth since I use Windows. Is R 3.6.2 the most recent update for Mac? I believe the most recent update is R 4.1.0, so maybe updating R + RStudio and then re-installing all your packages will fix it?

As far as updating, I believe there is some brms warnings in various vignettes and documentation about writing over an existing brms.fit object. So consider the following:

fit <- brm(Mean.HR ~ season + (1 | Date),
           data = turtledata.HRmean)

fit <- brm(Mean.HR ~ season + (1 | Date),
           data = turtledata.HRmean, iter = 4000)

In this case, Stan will attempt to compile the same model twice, and it will attempt to overwrite an object of the same name (i.e., fit). You can either run rm(fit) between the calls or do the second call like this (I think, the code may not be correct because I usually don’t use update()):

fit <- update(fit, iter = 4000)
7

I will definitely have to try updating R. I know that I do have the most updated R studio though.

I was doing some futher Googling and saw that I might need to download the app Xcode on Mac. It seems to take up a lot of space but I think that’s my next step, along with updating R also. Lots of updating to do!

Thank you again for your help!!

8

I installed the app Xcode and updated R and my computer to the newest version. This seemed to work and get rid of the sampling error.

I still have the divergent transitions error along with

Warning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low.
Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.

I used a simple model with mainly default settings:

m1double <- brm(double_model,
                data = turtledata.HRmean, seed = 12345)
9

Can you provide a plot of the outcome and perhaps give us some details of what the outcome is? I only found this:

The response variable is Mean.HR, season is a fixed effect, and Date and Turtle are random effects.

10

Here is a reproducible example of the data set:

>dput(head(turtledata.HRmean, 10))

structure(list(season = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L), .Label = c("beginning", "middle", "end"), class = "factor"), 
    Turtle = c("R3L1", "R3L1", "R3L1", "R3L1", "R3L1", "R3L1", 
    "R3L1", "R3L1", "R3L1", "R3L1"), Date = c("2015-05-24", "2015-05-25", 
    "2015-05-26", "2015-05-27", "2015-05-28", "2015-05-29", "2015-05-30", 
    "2015-05-31", "2015-06-01", "2015-06-02"), Cycle = 1:10, 
    Mean.HR = c(22.068574015748, 15.421802020202, 29.1950055335968, 
    28.754163099631, 22.1920775423729, 32.5265151832461, 26.0038786046512, 
    19.6938880769231, 10.6471096774194, 9.14997920792079), SD = c(16.8350435805238, 
    4.53600228857138, 16.38215204534, 14.8497129255194, 9.49196460499246, 
    17.1430413978331, 5.15496125729505, 4.53101262759239, 2.00900311400501, 
    1.59621326146081)), row.names = c(NA, -10L), groups = structure(list(
    season = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L), .Label = c("beginning", "middle", "end"), class = "factor"), 
    Turtle = c("R3L1", "R3L1", "R3L1", "R3L1", "R3L1", "R3L1", 
    "R3L1", "R3L1", "R3L1", "R3L1"), Date = c("2015-05-24", "2015-05-25", 
    "2015-05-26", "2015-05-27", "2015-05-28", "2015-05-29", "2015-05-30", 
    "2015-05-31", "2015-06-01", "2015-06-02"), .rows = structure(list(
        1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L), ptype = integer(0), class = c("vctrs_list_of", 
    "vctrs_vctr", "list"))), row.names = c(NA, -10L), class = c("tbl_df", 
"tbl", "data.frame"), .drop = TRUE), class = c("grouped_df", 
"tbl_df", "tbl", "data.frame"))

Basically, each Turtle has an average HR taken per day (this is what the Mean.HR value is). Each turtle was measured anywhere from 9-14 days. Here is the code for the double model that I used to run the brm function:

double_model = bf(Mean.HR ~ season + (1|Date), sigma ~ (1|Turtle))

Here is the output of the function brm:

summary(m1double)
 Family: gaussian 
  Links: mu = identity; sigma = log 
Formula: Mean.HR ~ season + (1 | Date) 
         sigma ~ (1 | Turtle)
   Data: turtledata.HRmean (Number of observations: 119) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Group-Level Effects: 
~Date (Number of levels: 78) 
              Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)     3.17      0.64     1.67     4.31 1.02
              Bulk_ESS Tail_ESS
sd(Intercept)      171      442

~Turtle (Number of levels: 11) 
                    Estimate Est.Error l-95% CI u-95% CI Rhat
sd(sigma_Intercept)     1.02      0.41     0.43     2.02 1.03
                    Bulk_ESS Tail_ESS
sd(sigma_Intercept)      125      160

Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept          16.51      1.51    13.15    19.15 1.01
sigma_Intercept     1.26      0.41     0.30     1.92 1.02
seasonmiddle       -0.08      1.72    -3.18     3.46 1.01
seasonend          -1.71      1.88    -5.15     2.17 1.01
                Bulk_ESS Tail_ESS
Intercept            530      589
sigma_Intercept      183      227
seasonmiddle         498      285
seasonend            524      806

I could use that information for my data analysis, but I’m just not sure if it’s correct because of all the error messages (which were listed above). Does this help?

11

Would you be so kind and plot Mean.HR?

12

of course. I’m not sure exactly what you are looking for (I’m just a beginner with this :) ). Would this work?

Screen Shot 2021-06-08 at 8.54.38 AM
Screen Shot 2021-06-08 at 8.54.38 AM1446×906 156 KB

13

So, one thing that strikes me is that the outcome can only be positive (is it a real number or is it a count?) Perhaps try a family which accounts for that?

14

And, we love beginners :)

15

The Mean.HR is a real number like 14.56 or 17.8. Things like that. What do you mean by family?

16

The family = ... argument in brms allows you to specify different kinds of likelihood that could describe your data. In your case, you’re estimating a Gaussian (normal) likelihood or family, which implies that your data can be positive or negative. There are other families that you can estimate that may make more sense with your data. For example, both the family = Gamma and family = lognormal specifications are likelihoods that assume strictly positive continuous variables.

The brms documentation for families is really good, and you can find it here. There is also a lot of really good information about the various types of likelihood families that Stan generally supports. The current Stan Functions Reference guide groups these by the kinds of families, and this link will take you to the section on positive continuous distributions.

-----Edit (Additional Context)-----
The reason that potentially changing your family argument might help in your case is that there are parts of the posterior that do no need to be explored. The issue with Bayesian estimation is that we aren’t guaranteed a nice, well-defined posterior, meaning that sometimes our tools (like Hamiltonian Monte Carlo in the case of Stan) will fail to give us good descriptions of the posterior distribution where we make our inferences. We can help avoid that by doing things that “control” the posterior distribution. Historically, this was done with the use of conjugate priors, but now we have other more useful/flexible tools. Since you can’t observe a negative value of average heart rate, specifying a likelihood that only allows positive values ensures that we don’t try to estimate any values that don’t make sense. Similarly, using more regularizing priors can help to keep the posterior well-behaved.

1 Like