Heritability from betabinomial GLMM variance estimates using pedigree data

Interfaces brms
1

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 (h2) 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, VP) is composed of the genetic variance (VG) and environmental variance (VE), so that
VP = VG + VE

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

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

2 Likes
2

See if the discussion here helps you.

3

Both that discussion and this one are related to my questions. Phylogenetic models are analogous to pedigree modes, and upon further reading I think the best way to address my second question will be to use the hypothesis framework from the phylogenetic model vignette.

I am still trying to figure out whether to use the identity phi to calculate residual variance (𝛔2𝛜) when the genetic variance (𝛔2𝛂) is estimated on the logit link scale.

4

The phi in the output is the actual phi as you don’t predict it (but assume it to be constant across observations) and hence brms uses no link function (i.e., an identity link).

5

So from my betabinomial models I should be able to calculate residual variance directly as 𝛟𝛑2/3?

6

If that is the right variance formula then yes, but i haven’t checked its correctness.