[R-sig-ME] Random effect variance estimates not matching in MCMCglmm and brms

[Simon Tapper]

Hi All,

I've run an analysis in MCMCglmm and in brms and am getting somewhat
different results in the random slope and sigma estimates between the two
packages. For my analysis, I applied a treatment to a group of animals, and
am trying to estimate heterogeneous variance components for each treatment.
In both packages, I've run the model with uninformative priors for the
random effect and residual components. I'm wondering if the discrepancy in
results is from how the two packages estimate the random slopes, or if from
differences in how I set the priors. I cannot provide reproducible data for
this question, but perhaps one can help me directly translate the MCMCglmm
priors into brms.

MCMCglmm model:

prior_MCMC <- list(B = list(mu = c(35, 0, -0.5, 0.05, 0), V=diag(5) * 1.5),
                  R = list(V = diag(2), nu = 1.002),
                  G = list(G1 = list(V = diag(2), nu =1.002),
                              G2 = list(V = diag(2), nu = 1.002)))

parallel::mclapply(1:3, function(i) {
  MCMCglmm(DV ~ Treatment*X + Y,
                 random =~us(at.level(Treatment, 1)+
                               at.level(Treatment, 1):X):ID +
                 us(at.level(Treatment, 2)+
                      at.level(Treatment, 2):X):ID,
                 data=data.frame(data),
                 rcov = ~idh(Treatment):units,
                 prior=prior_MCMC,
                 pr=T,
                 nitt=120000,
                 burnin=40000,
                 thin=30,
                 verbose=F)
}, mc.cores=4)

From my understanding, in MCMCglmm, when V=1, this is equivalent to an
inverse gamma distribution with shape and scale parameters set to nu/2. So,
whe nu = 1.002 in MCMCglmm, this would be 0.501 for a brms inverse gamma
prior. For the fixed effects in MCMCglmm, are the priors set properly for
mu? where DV = 35, Treatment_1 = 0, Treatment_2 = -0.5, X= 0.05, Y = 0.

brms model:

prior_treat <- c(
  set_prior("normal(35, 1.5)", class = "Intercept"),
  set_prior("normal(-0.5, 1.5)", class="b", coef="Treatment"),
  set_prior("normal(0.05, 1.5)", class="b", coef="X"),
  set_prior("normal(0, 1.5)", class="b", coef="Y"),
  set_prior("inv_gamma(0.501, 0.501)", class = "sd"),
  set_prior("inv_gamma(0.501, 0.501)", class = "b", lb=0, dpar="sigma"))

brm(bf(DV ~ Treatment*X + Y +
               (1 + X|gr(ID, by=Treatment)),
                 sigma~Treatment),
                data=data.frame(data),
                warmup = 1000,
                iter = 10000,
                chains = 3,
                thin = 10,
                cores = 4,
                seed = 123,
                sample_prior = TRUE,
                prior = prior_treat,
                control = list(adapt_delta = 0.95, max_treedepth=13))

Results:
In MCMCglmm,the variance for the slope estimates (for both treatments) ±
95% credible intervals are ~0.2 (0.05-0.5) and sigma = ~0.1. In brms, the
variance for the slope estimates (for both treatments) ± 95% credible
intervals are one order of magnitude lower, ~0.01 (0.0025-0.01) and sigma =
0.01. Oddly, the variance and credible interval estimates for the random
intercepts appear pretty similar between packages.

Any help would be appreciated,
Simon

	[[alternative HTML version deleted]]