Different prior settings for within and between-group variance

Interfaces brms
, , ,
1

Hi, I was just interested to know if anyone had any insight into a model I am currently trying to run. In brief I am interested in seeing how the mean and variance in behaviour (in this case animals) might be correlated to reproductive success. In my example I have multiple measurements of an individuals behaviour within a given year but only a singular measure of their reproductive success.

A mock dataset for what I am looking at is below.
|Column 1 | Column 2 | Column 3 | Column 4 | E | F|
|Ind_ID | Foraging_time | Temperature | Year | Breed_success | Sub|
|Animal_1 | 40 | 5 | 2023 | 1 | T|
|Animal_1 | 20 | 5 | 2023 | NA | F|
|Animal_1 | 30 | 7 | 2023 | NA | F|
|Animal_1 | 25 | 7 | 2023 | NA | F|
|Animal_2 | 10 | 5 | 2023 | 0 | T|
|Animal_2 | 15 | 5 | 2023 | NA | F|
|Animal_2 | 10 | 7 | 2023 | NA | F|
|Animal_2 | 15 | 7 | 2023 | NA | F|
|Animal_3 | 60 | 5 | 2023 | 1 | T|
|Animal_3 | 30 | 5 | 2023 | NA | F|
|Animal_3 | 20 | 7 | 2023 | NA | F|
|Animal_3 | 40 | 7 | 2023 | NA | F|

I am attempting to run this model in brms with mean and variance components for behaviour and just mean for the reproductive success. I would like to run this as a bivariate model so that I can account for “fixed” effects on both response variables. What I am wondering is if it is possible to set a prior so that the within-individual variance of a random effect can be fixed at 0 whilst still allowing the between-individual variance to be calculated in order to quantify the correlation between foraging and reproductive success. I have attempted to use the constant(0.0001) setting when setting the prior for the breed_success model, but I think it is removing all variance and therefore not allowing the calculation of the correlation between the two measures (when the model runs, the variance is set to 0 but the correlations are also remarkably close to 0 for me). I have previously run a similar model in MCMCglmm before that allows the within- and among-group variances to be set differently (hence the within- can be fixed to 0) but do not know if there is a similar way of doing this in brms. An example of the type of prior I am thinking of in MCMCglmm can be found on PAGE 17 of this document https://tomhouslay.files.wordpress.com/2017/02/indivvar_plasticity_tutorial_mcmcglmm1.pdf. If anyone knows if this is possible for a prior setting to allow there to be among-individual variance without any within-individual that would be greatly appreciated.

The model structure currently would look like this

m.f.t ← bf(Foraging_time ~ Temperature + Year + (1|q|Ind_ID)) + gaussian()
m.b.s ← bf(Breed_success | subset(Sub) ~ Year + (1|q|Ind_ID)) + bernoulli()

model ← brm(m.f.t + m.b.s,
data = data,
prior = prior,
warmup = 1000,
iter = 5000,
thin=4,
chains = 4,
seed = 12345)

Although, I understand I only have single measures of an individuals Reproductive success and therefore setting a random effect is odd, it is required to calculate the correlation between the two responses.

Any help on this topic would be greatly appreciated and I look forward to hearing from anyone!
Cheers,
Freddie

2

I’m not an expert so take this with a grain of salt, but I am not sure this approach works with a Bernoulli outcome. The “within-group” variation is the residual variation in your model. To mimic the MCMCglmm tutorial you would set a prior on sigma, the residual standard deviation in brms. However, the tutorial has a Gaussian outcome (fitness), while you have a Bernoulli outcome (breeding success), and no residual variance is estimated. So I don’t think you can do this in MCMCglmm either?

3

Hi and thanks for replying! That does make sense to me, hadn’t really thought about the fact with their being no estimated residual variance with a bernoulli variance. If I were therefore able to have a different breeding success variable that would follow a gaussian distribution, on what level of sigma would I set the prior on?
Would I have to include a residual variance part of the model as in
b.s.gaussian ← bf(gauss.fitness ~ Year + (1|Ind_ID), sigma ~ (1|Ind_Id) + gaussian() and then set a constant prior on the Ind_ID sigma component of the fitness response?
(set_prior(“constant(0.0001)”, class = “sd”, group = “Ind_Id”, resp = “gauss.fitness”, dpar = “sigma”)) ?
I guess this would solely limit the within-individual variance in to 0?

The other option I guess would be having bf(gauss.fitness ~ Year + (1|Ind_ID)) + gaussian() and then having a prior set_prior(“constant(0.0001)”, class = “sigma”, resp = “gauss.fitness”). I guess this would mean setting all variance to 0 rather than just the within individual variance due to individual ID?

If you could offer any insight that would be amazing!

4

I think your second option. I don’t think modelling sigma ~ (1 | Ind_Id) makes much sense, sigma is already the standard deviation among individuals. So the equivalent model to the MCMCglmm tutorial would be (I think):

m.f.t <- bf(Foraging_time ~ Temperature + Year + (1|q|Ind_ID)) + gaussian()
m.b.s <- bf(Fitness | subset(Sub) ~ Year + (1|q|Ind_ID)) + gaussian()

With a small constant prior on sigma for Fitness to move all the variance into the random effect. I would recommend trying the approach on the data from the tutorial to see if it actually works the way we think it does.

Not sure if I follow. That prior just sets the residual standard deviation to a small value, it does not affect the variance captured by the random effect.

5

Hi again,
So you are right with this. I have been playing around with the data from the tutorial and the model fitting a small constant prior on sigma for fitness produces similar parameter estimates for the covariances and correlations.

MCMCglmm model:

prior_biv_RR_px <- list(R = list(V = diag(c(1,0.0001),2,2), nu = 0.002, fix = 2),
G = list(G1 = list(V = matrix(c(1,0,0,0,1,0,0,0,1),3,3,byrow = TRUE),nu = 3,alpha.mu = rep(0,3),alpha.V = diag(25^2,3,3))))

mcmc_biv_RR <- MCMCglmm(cbind(scale(aggression),
fitness) ~ trait-1 +
at.level(trait,1):opp_size +
at.level(trait,1):scale(assay_rep, scale=FALSE) +
at.level(trait,1):block +
at.level(trait,1):scale(body_size),
random =~ us(trait + opp_size:at.level(trait,1)):ID,
rcov =~ idh(trait):units,
family = c("gaussian","gaussian"),
prior = prior_biv_RR_px,
nitt=950000,
burnin=50000,
thin=450,
verbose = TRUE,
data = as.data.frame(df_plast),
pr = TRUE,
saveX = TRUE, saveZ = TRUE)

brms.model:

agg.mod <- bf(scale(aggression) ~ opp_size + 
                scale(assay_rep, scale = FALSE) +
                block +
                scale(body_size) +
                (1+opp_size|q|ID)) + gaussian()

fit.mod <- bf(rel.fitness|subset(fitness.sub) ~ (1|q|ID)) + gaussian()

brm_biv_prior <- c(set_prior("lkj(1)", class = "cor"),
                   set_prior(student_t(3, 1, 2.5), class = "Intercept", resp = "relfitness"),
                   set_prior("student_t(3, 0, 2.5)", class = "sd", resp = "relfitness", group = "ID"),
                   set_prior("constant(0.0001)", class = "sigma", resp = "relfitness"),
                   set_prior("(flat)", class = "b", resp = "scaleaggression"),
                   set_prior("student_t(3, 0, 2.5)", class = "Intercept", resp = "scaleaggression"),
                   set_prior("student_t(3,0,2.5)", class = "sd", resp = "scaleaggression", group = "ID"),
                   set_prior("student_t(3,0,2.5)", class = "sigma", resp = "scaleaggression")
                   )

brm_biv_RR <- brm(agg.mod + fit.mod + set_rescor(FALSE),
                   data = df_plast,
                   prior = brm_biv_prior,
                   warmup = 1000,
                   iter = 3000,
                   thin = 2,
                   chains = 4,
                   init = "random",
                   seed = 12345)

I have kept all the default priors for the brms model apart from the sigma for now.

Although this model produces estimates similar to what we would expect based on the MCMCglmm model, I have to take these pretty cautiously as the model does not mix that well at all. there are a lot of transitions that exceed max tree depth and the sample sizes produced are quite low.

Think I need to play around with it a bit more as it could be a sample size issue with the fewer points for fitness but if you have any suggestions of course I will take them on board!
Thanks for all your help so far!

6

Whoops, turns out just more iterations and thinning was the key…
Thanks again @Ax3man . Now I know it works for the tutorial dataset just need to run it on my own data again!

1 Like