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

[Jarrod Hadfield]

Hi,

With a 2x2 us structure the marginal dsitribution of a variance is
inverse-Wishart with nu*=nu-1 and V*=V*nu/(nu-1), so in your case
nu*=0.002 and V*=501. This is inverse gamma with shape = nu*/2= 0.001
and scale (nu*)*(V*)/2=0.501, not 0.501 for both.

With an idh structure the marginal dsitribution of a variance is simply
inverse-Wishart with nu*=nu and V*=V which is inverse-gamma with a shape
and scale (in your case) of 0.001.

I don't use brms so I'm not exactly clear what model/prior is being
used, but is the second argument to normal in the prior a variance?
Typically a standard deviation is used so you might want to square-root
the 1.5 in the fixed effect prior.

I'm also not sure what is happening about the intercpet-slope covariance
in brms?

Cheers,

Jarrod



On 01/04/2021 20:48, Simon Tapper wrote:
> This email was sent to you by someone outside the University.
> You should only click on links or attachments if you are certain that the email is genuine and the content is safe.
>
> 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]]
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336. Is e buidheann carthannais a th’ ann an Oilthigh Dhùn Èideann, clàraichte an Alba, àireamh clàraidh SC005336.