Is it possible to get an error with rhat() due to convergence problems?

Interfaces brms
,
1

As a general context, within my doctoral thesis I am doing a simulation study to evaluate the statistical characteristics of the Bayesian and frequentist framework in the application of multilevel models to data from single-case designs. This involves simulating data with four models varying the number of random effects. We are currently trying to replicate one condition within the 144 we have because we noticed that many missing values appear. Specifically, this condition tries to simulate data with a straightforward model that has no random effects and these data are analyzed with four models that vary in random effects and then compare the performance of different fit indices (frequentist and Bayesian). In addition, this condition has a large effect size (d = 2.70) and 10 repeated measures per participant, considering that there are 5 per simulated dataset.

We have 500 replicates per condition. When replicating this condition I see that the following problem appears when I want to evaluate the convergence of the model:

Error in UseMethod(“rhat”) :
no applicable method for ‘rhat’ applied to an object of class “brmsfit”.

I was wondering if this error can come up because the model has not converged and then there is no sense in calculating a rhat. Attached is the code of the model and the fitting result. If you need any additional info, I can also try to provide it. Thank you very much in advance.

fit_mxc2 <- brm(y ~ x + (1+x|id), data = databrms,
                        prior = c(set_prior("normal(0,1000000)", class = "b"),
                                  set_prior("cauchy(0,20)", class = "sd"),
                                  set_prior("lkj(2)", class="cor"),
                                  set_prior("inv_gamma(0.001,0.001)", class = "sigma")),
                        warmup = 400,
                        iter = 1000, chains = 2, control = list(adapt_delta = 0.95), cores = 4)

fit_max_c2 = update(fit_mxc2,newdata=databrms, cores = 4)

The results I have obtained are the following:

fit_max_c2
Family: gaussian
Links: mu = identity; sigma = identity
Formula: y ~ x + (1 + x | id)
Data: databrms (Number of observations: 50)
Draws: 2 chains, each with iter = 1000; warmup = 400; thin = 1;
total post-warmup draws = 1200

Multilevel Hyperparameters:
~id (Number of levels: 5)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.42 0.40 0.02 1.38 1.00 356 586
sd(x) 0.53 0.49 0.01 1.80 1.01 306 372
cor(Intercept,x) -0.08 0.45 -0.88 0.77 1.00 862 657

Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 4.82 0.35 4.06 5.57 1.00 738 589
x 2.92 0.45 2.03 3.77 1.00 806 490

Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 1.16 0.12 0.94 1.43 1.00 1639 729

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).

However, if I want to use rhat(), the following error appears:

rhat(fit_max_c2)
Error in UseMethod(“rhat”) :
no applicable method for ‘rhat’ applied to an object of class “brmsfit”

  • Operating System: Windows 11 Home
  • brms Version: 2.21.0

Thanks again.

2

There is currently no method called “rhat” for brmsfit objects. If you want to compute rhat, you may want to look at, for example, the summarize_draws function of the posterior package. E.g.

draws ← as.draws.array(brmsfit)
posterior::summarize_draws(draws)

3

Thank you very much, Paul, for your prompt reply! With the function you are telling me about, I think I will be able to do what I need to do. However, I wanted to tell you that when testing it on much simpler models (reproducible example step) the rhat() function does return a result. Have you made this change recently? :)

# Data simulation 
set.seed(123)  # For reproducibility matters
n <- 100  # Number of observations
x <- rnorm(n, mean = 5, sd = 2)  # Independent variable
y <- 3 + 2 * x + rnorm(n, mean = 0, sd = 1)  # Dependent variable
my_data <- data.frame(x = x, y = y)

# Fitting the model with brms
fit <- brm(formula = y ~ x, data = my_data, family = gaussian(), chains = 2, iter = 2000)

# Summary of results
summary(fit)

# Extract r hat values
rhat_values <- rhat(fit)
print(rhat_values)

summary(fit)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: y ~ x
Data: my_data (Number of observations: 100)
Draws: 2 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 2000

Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 3.04 0.30 2.46 3.62 1.00 1552 1197
x 1.97 0.05 1.87 2.08 1.00 1509 1276

Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.98 0.07 0.85 1.13 1.00 2182 1478

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).

rhat_values ← rhat(fit)
print(rhat_values)
b_Intercept b_x sigma Intercept lprior lp__
1.000349 1.000360 1.003052 1.001198 1.001469 1.000092

4

There is an rhat extractor functon in the bayesplot package. Are you maybe using that one?