*I’m a bit new to brms and bayesian modeling in general so this might be an obvious one.*

I am trying to calculate heritability (h^{2}) of a single trait using a pedigree. In a simplified framework, I’m interested in what proportion of the variance of the trait of interest (phenotypic variance, V_{P}) is composed of the genetic variance (V_{G}) and environmental variance (V_{E}), so that

V_{P} = V_{G} + V_{E}

Heritability is calcuated as the proportion of VP attributable to additive or sex-linked genetic variance.

h^{2} = V_{G} /V_{P}

This calculation is analogous to that for repeatability (R ) à la Nakagawa and Schielzeth (2010).

R = 𝛔^{2}_{𝛂} / (𝛔^{2}_{𝛂} + 𝛔^{2}_{𝛜})

where 𝛔^{2}_{𝛂} is the between group variance and 𝛔^{2}_{𝛜} is the residual variance.

I fitted a beta binomial GLMM (following the Custom Response Distributions vignette), where Amatrix97 is the genetic relatedness matrix of the random effect IDA.

`myModel_betaBinomial <- brm(formula = trt_eggs|vint(total_eggs) ~ Gen + (1|gr(IDA, cov = Amatrix97)), data=PM1997_brms, family = beta_binomial2, stanvars = stanvars, data2 = list(Amatrix97 = Amatrix97), warmup = 1000, iter = 20000, chains = 1, inits= "0", control = list(max_treedepth = 20, adapt_delta = 0.95))`

Question 1.

The residual variance of the beta binomial distribution is the product of the overdispersion parameter (𝛟) and the distribution-specific variance:

𝛔^{2}_{𝛜} = 𝛟𝛑^{2}/3

The brms output states that while 𝛍 is estimated with the logit link, 𝛟 is estimated with the identity link. Can I use 𝛟 as estimated by the model to calculate the residual variance?

Question 2.

Assuming the answer to question 1 is yes, I would calculate variance components and heritability in R as:

`dsv <- (pi^2)/3 #distribution-specific variance`

`posterior_mybb <- posterior_samples(myModel_betaBinomial, pars = c("sd_IDA", "phi")) %>% mutate(Va = sd_IDA__Intercept^2, Ve = phi*dsv) %>% select(Va, Ve)`

`h2a_mybb <- posterior_mybb[,1] /(posterior_mybb[,1]+ posterior_mybb[,2])`

`est_h2a_mybb <- mean(h2a_mybb)`

What is the best approach for calculating an interval (like HPDinterval in MCMCglmm) around this estimate?

Thanks in advance!

R version 4.0.2 (2020-06-22)

Platform: x86_64-apple-darwin17.0 (64-bit)

Running under: macOS Catalina 10.15.4

Brms 2.13.3